Commit 00bb10d2 authored by Florent Chatelain's avatar Florent Chatelain
Browse files

fix typos

parent e76dcb05
......@@ -3,43 +3,44 @@
This notebook can be run on mybinder: [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/git/https%3A%2F%2Fgricad-gitlab.univ-grenoble-alpes.fr%2Fchatelaf%2Fml-sicom3a/master?urlpath=lab/tree/notebooks/7_Clustering/N1_Kmeans_basic.ipynb)
%% Cell type:markdown id: tags:
# KMEANS basics
The purpose of this lab is to implement simple 1D Kmeans clustering algorithm, and compare the obtained results with those obtained using skleran implementation
The purpose of this lab is to implement simple 1D Kmeans clustering algorithm, and compare the obtained results with those obtained using sklearn implementation
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
%matplotlib inline
```
%% Cell type:markdown id: tags:
## import data from matlab file :
%% Cell type:code id: tags:
``` python
Data=loadmat('fictitious_train.mat')
Data = loadmat("fictitious_train.mat")
print(Data.keys())
X=Data.get('Xtrain')
print('dim of X:{}'.format(X.shape))
X = Data.get("Xtrain")
print("dim of X:{}".format(X.shape))
```
%% Cell type:markdown id: tags:
## Compute the histogram
%% Cell type:code id: tags:
``` python
bins=np.arange(np.min(X)-1,np.max(X)+2,1)
hist_val,bins=np.histogram(X, bins=bins)
bins = np.arange(np.min(X) - 1, np.max(X) + 2, 1)
hist_val, bins = np.histogram(X, bins=bins)
print(hist_val)
print(bins)
```
%% Cell type:markdown id: tags:
......@@ -47,13 +48,13 @@
### or directly visualize the histogram
%% Cell type:code id: tags:
``` python
bins=np.arange(np.min(X)-1,np.max(X)+2,1)
plt.scatter(X,np.zeros_like(X)+.5,c='red',marker='+')
n,bin_edges,patches=plt.hist(x=X,bins=bins, color='blue',histtype='step')
bins = np.arange(np.min(X) - 1, np.max(X) + 2, 1)
plt.scatter(X, np.zeros_like(X) + 0.5, c="red", marker="+")
n, bin_edges, patches = plt.hist(x=X, bins=bins, color="blue", histtype="step")
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWsAAAD4CAYAAAAqw8chAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAALWElEQVR4nO3dYYhl91nH8d/jbqRtWmkhUTSb3akg1VJoIkOsBiSmRaIt9W0K7QsR9k3VVArF+s73UuoLEZa0KjS2SJqCBK0NtJtS0OhsEjXpplDibhsTyQRpmvjCmvbxxcxmZ3b37tzZuXfO/Hc+Hxgys/fcc56cmfvds2fP2VvdHQAOth+begAAdibWAAMQa4ABiDXAAMQaYABHl7HSm266qVdWVpaxaoDr0pkzZ17q7ptnPb6UWK+srGRtbW0Zqwa4LlXV+as97jQIwADEGmAAYg0wALEGGIBYAwxArAEGMNele1V1LskrSX6Y5LXuXl3mUABst5vrrH+tu19a2iQAzOQ0CMAA5o11J/lKVZ2pqpNXWqCqTlbVWlWtra+vL25CJreyklRN/+FfMOAwm/c0yJ3d/XxV/WSSR6rqme7++tYFuvtUklNJsrq66u1nriPnzycH4Q2FqqaeAKYz15F1dz+/+d8Xk3wpyR3LHAqA7XaMdVXdWFVvufB5kl9P8tSyBwPgonlOg/xUki/Vxp9Bjyb56+7+8lKnAmCbHWPd3c8mefc+zALADC7dAxiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQYg1gADEGuAAYg1wADEGmAAYg0wALEGGIBYAwxArAEGINYAAxBrgAGINcAAxBpgAGINMACxBhiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTCAuWNdVUeq6omqeniZAwFwud0cWd+X5OyyBgFgtrliXVXHkrw/yf3LHQeAK5n3yPrTST6R5EezFqiqk1W1VlVr6+vrCxluSisrSdW0HysrU+8FruQg/Gz4+Th8ju60QFV9IMmL3X2mqu6atVx3n0pyKklWV1d7YRNO5Pz5pCf+v6iadvtc2UH42Uj8fBw28xxZ35nkg1V1LskXktxdVZ9b6lQAbLNjrLv7k919rLtXktyb5Kvd/eGlTwbA61xnDTCAHc9Zb9Xdp5OcXsokAMzkyBpgAGINMACxBhiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQYg1gADEGuAAYg1wADEGmAAYg0wALEGGIBYAwxArAEGINYAAxBrgAGINcAAxBpgAGINMACxBhiAWAMMQKwBBiDWAAMQa4ABiDXAAHaMdVW9oar+uar+taqerqo/3o/BALjo6BzL/G+Su7v71aq6Ick3qurvu/ufljwbAJt2jHV3d5JXN7+8YfOjlzkUANvNdc66qo5U1ZNJXkzySHc/ttyxANhqrlh39w+7+7Ykx5LcUVXvunSZqjpZVWtVtba+vr7oOYFLnDiRVE3/sbIy9Z44HHZ1NUh3fy/J6ST3XOGxU9292t2rN99884LGA2Y5dy7pnv7j/Pmp98ThMM/VIDdX1Vs3P39jkvcleWbZgwFw0TxXg/x0kr+qqiPZiPvfdPfDyx0LgK3muRrk35Lcvg+zADCDOxgBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQYg1gADEGuAAYg1wADEGmAAYg0wALEGGIBYAwxArAEGINYAAxBrgAGINcAAxBpgAGINMACxBhiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQawY6yr6taq+lpVna2qp6vqvv0YDICLjs6xzGtJPt7dj1fVW5KcqapHuvubS54NgE07Hll39wvd/fjm568kOZvklmUPBsBF8xxZv66qVpLcnuSxKzx2MsnJJDl+/PgCRuPEiaRq6ik25gCmNXesq+rNSb6Y5GPd/f1LH+/uU0lOJcnq6movbMJD7Ny5qScADoq5rgapqhuyEeoHuvuh5Y4EwKXmuRqkknwmydnu/tTyRwLgUvMcWd+Z5CNJ7q6qJzc/fnPJcwGwxY7nrLv7G0kOwF9zARxe7mAEGIBYAwxArAEGINYAAxBrgAGINcAAxBpgAGINMACxBhiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQYg1gADEGuAAYg1wADEGmAAYg0wALEGGIBYAwxArAEGINYAAxBrgAGINcAAxBpgAGINMIAdY11Vn62qF6vqqf0YCIDLzXNk/ZdJ7lnyHABcxY6x7u6vJ/nvfZgFgBmOLmpFVXUyyckkOX78+KJWC687cSKpmnqKjTm46CB8X06cSM6dm3aGZVtYrLv7VJJTSbK6utqLWi9ccL2/GEd1EL4vU/9msR9cDQIwALEGGMA8l+59Psk/JnlHVT1XVb+z/LEA2GrHc9bd/aH9GASA2ZwGARiAWAMMQKwBBiDWAAMQa4ABiDXAAMQaYABiDTAAsQYYgFgDDECsAQYg1gADWNibD1xX7roryek9PHeL03Ou58Lz5l1+UdvdzfqffDK57bbLH1vktva6H/Z7e9fy/FnPmXddu9nmbuebZ/llfo/2+/s/EEfWAANwZL3Vhd/VH310+9e7OYK58NxLf33WOvayzb1sd14XjqhffvnK21nUtva6H/Z7e9fy/FnPuWCnde1mm7udb57ll/k92u/v/4AcWQMMoLoX/962q6urvba2tvD17pu77ko9ejrXtGucs977dha9zmVuzznrxbrGdVfl2l6vB0hVnenu1VmPO7IGGIAj6xmuh9+p4bC4Hl6vjqwBrgNiDTAAsQYYgFgDDECsAQYg1gADEGuAAYg1wADEGmAAYg0wALEGGIBYAwxArAEGINYAA5gr1lV1T1V9q6q+XVV/uOyhANhux1hX1ZEkf5bkN5K8M8mHquqdyx4MgIvmObK+I8m3u/vZ7v5Bki8k+a3ljgXAVvO8u/ktSb675evnkvzSpQtV1ckkJze/fLWqvrX38aZVteun3JTkpcVPMiT7Yjv7Y7uF749reL0eFBf2xYmrLTRPrK+0Cy57A53uPpXk1FyjXaeqau1qb8tzmNgX29kf29kfF827L+Y5DfJcklu3fH0syfPXOhgAuzdPrP8lyc9V1dur6seT3Jvkb5c7FgBb7XgapLtfq6rfTfIPSY4k+Wx3P730ycZ0qE8DXcK+2M7+2M7+uGiufVE9+vu3AxwC7mAEGIBYAwxArBfA7fgXVdWtVfW1qjpbVU9X1X1TzzS1qjpSVU9U1cNTzzK1qnprVT1YVc9s/oz88tQzTamq/mDzdfJUVX2+qt4wa1mx3iO341/mtSQf7+5fSPKeJB895PsjSe5LcnbqIQ6IP03y5e7++STvziHeL1V1S5LfT7La3e/KxgUc985aXqz3zu34W3T3C939+Obnr2TjxXjLtFNNp6qOJXl/kvunnmVqVfUTSX41yWeSpLt/0N3fm3aqyR1N8saqOprkTbnKPSxivXdXuh3/0MZpq6paSXJ7ksemnWRSn07yiSQ/mnqQA+Bnk6wn+YvN00L3V9WNUw81le7+zyR/kuQ7SV5I8nJ3f2XW8mK9d3Pdjn/YVNWbk3wxyce6+/tTzzOFqvpAkhe7+8zUsxwQR5P8YpI/7+7bk/xPkkP7dzxV9bZs/Cn87Ul+JsmNVfXhWcuL9d65Hf8SVXVDNkL9QHc/NPU8E7ozyQer6lw2To/dXVWfm3akST2X5LnuvvAnrQezEe/D6n1J/qO717v7/5I8lORXZi0s1nvndvwtqqqycU7ybHd/aup5ptTdn+zuY929ko2fi69298wjp+tdd/9Xku9W1Ts2f+m9Sb454UhT+06S91TVmzZfN+/NVf7CdZ5/dY+rcDv+Ze5M8pEk/15VT27+2h91999NOBMHx+8leWDzwObZJL898TyT6e7HqurBJI9n4yqqJ3KVW8/dbg4wAKdBAAYg1gADEGuAAYg1wADEGmAAYg0wALEGGMD/A/6q86ZPadgCAAAAAElFTkSuQmCC)
......@@ -61,92 +62,91 @@
%% Cell type:markdown id: tags:
## Implementation of Kmean on a simple case
In this example, the number of clusters is assumed to be known.
### Exercize 1 :
### Exercise 1 :
- Explain/ comment the code below
- What is the main problem left aside by this code?
%% Cell type:code id: tags:
``` python
K=2 #nb of clusters
p=1 # dimension (the code below is given for p=1 only)
K = 2 # nb of clusters
p = 1 # dimension (the code below is given for p=1 only)
```
%% Cell type:code id: tags:
``` python
N=X.size
idx=np.zeros((N,1))
muvec=np.zeros((K,1))
N = X.size
idx = np.zeros((N, 1))
muvec = np.zeros((K, 1))
change = True # Defines the test variable for the loop.
# Default is true (meaning that a new iteration will be performed
perm = np.random.permutation(N)[0:2]
# takes two different random integers between 0 and $N$
change = True # Defines the test variable for the loop.
# Default is true (meaning that a new iteration will be performed
perm=np.random.permutation(N)[0:2]
# takes two different random integers between 0 and $N$
for k in range(0, K):
muvec[k] = X[perm[k], :] # Initialization of the cluster representatives (centers)
for k in range (0,K):
muvec[k] = X[perm[k],:] #Initialization of the cluster representatives (centers)
for i in range(0, N):
d = (X[i] - muvec) ** 2 # Computation of distances wrt cluster centers
idx[i] = np.where(d == d.min())[0] # label = index of closest center
for i in range (0,N):
d=(X[i] - muvec )**2 #Computation of distances wrt cluster centers
idx[i]=np.where(d==d.min())[0] #label = index of closest center
while change:
change=False
#update
for k in range (0,K):
muvec[k]= np.mean( X[idx == k] ) #compute new centers
#prediction
for i in range (0,N):
d=(X[i] - muvec )**2 #Computation of distances wrt cluster centers
index=np.where(d==d.min())[0]##label = index of closest center
if index != idx[i]: #check if some indices changed
change=True
idx[i]=index #replaces new index set
X0=X[idx==0]
X1=X[idx==1]
bins=np.arange(np.min(X)-1,np.max(X)+2,1)
n,bin_edges,patches=plt.hist(x=X,bins=bins, color='blue',histtype='step')
plt.scatter(X0,np.zeros_like(X0)+.5,c='red',marker='+', label='class 0')
plt.scatter(X1,np.zeros_like(X1)+.5,c='green',marker='+',label='class 1')
change = False
# update
for k in range(0, K):
muvec[k] = np.mean(X[idx == k]) # compute new centers
# prediction
for i in range(0, N):
d = (X[i] - muvec) ** 2 # Computation of distances wrt cluster centers
index = np.where(d == d.min())[0] ##label = index of closest center
if index != idx[i]: # check if some indices changed
change = True
idx[i] = index # replaces new index set
X0 = X[idx == 0]
X1 = X[idx == 1]
bins = np.arange(np.min(X) - 1, np.max(X) + 2, 1)
n, bin_edges, patches = plt.hist(x=X, bins=bins, color="blue", histtype="step")
plt.scatter(X0, np.zeros_like(X0) + 0.5, c="red", marker="+", label="class 0")
plt.scatter(X1, np.zeros_like(X1) + 0.5, c="green", marker="+", label="class 1")
plt.legend()
h=plt.gcf()
h = plt.gcf()
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWsAAAD4CAYAAAAqw8chAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQBUlEQVR4nO3df2xd5X3H8c+3iSdTyNQRRyPkJnYQqFCc2KlMQhYp8eKykLWACAhRKJGniUioIylUQesSQRaIxB+o6iQmIlPaLCJq+VGGgG3QSqlTgQbMBhpIQ+WSJfVNM2EcJSQLGYF+98e9Tvzbx/Y9Pvdrv1/SVXx9n/Oc731y78ePH59zrrm7AADl7QtZFwAAGBlhDQABENYAEABhDQABENYAEMD0NDqtqqrympqaNLoGgEmpvb39I3efNdTjqYR1TU2N2tra0ugaACYlMzs03OMsgwBAAIQ1AARAWANAAKmsWQOAJJ05c0b5fF6nT5/OupSyUVlZqVwup4qKilFtR1gDSE0+n9eMGTNUU1MjM8u6nMy5u7q7u5XP5zV//vxRbcsyCIDUnD59WjNnziSoi8xMM2fOHNNvGolm1mZ2UNIJSZ9L+szdG0a9JwBTEkHd11jHYzTLIH/p7h+NaS8AgHFhGQTAlLNlyxY98sgjqfTd3t6uBQsW6NJLL9X69etVqs8MSBrWLunnZtZuZusGa2Bm68yszczaurq6SlIcykNNjWSW/Y0rGCCCu+66Sy0tLero6FBHR4defvnlkvSbNKyXuftXJa2W9G0zW96/gbu3uHuDuzfMmjXk6e0I6NAhyT3726FhT8bFpNHYWLiVyM6dO7Vw4ULV1dXpjjvuGPD4448/rquuukp1dXW66aabdOrUKUnSM888o9raWtXV1Wn58kLk7du3T4sXL1Z9fb0WLlyojo6OPn0dOXJEH3/8sZYuXSoz09q1a/X888+X5HkkWrN29z8U//3QzP5V0mJJvypJBQCQkn379mnbtm167bXXVFVVpaNHjw5os2bNGt15552SpM2bN+uJJ57Q3Xffra1bt+qVV17RnDlzdOzYMUnS9u3btWHDBt1+++369NNP9fnnn/fp6/Dhw8rlcmfv53I5HT58uCTPZcSwNrPzJX3B3U8Uv/4rSVtLsncA6NEzm96zp+/91tYxd7l7927dfPPNqqqqkiRdeOGFA9q899572rx5s44dO6aTJ09q1apVkqRly5apublZt9xyi9asWSNJWrp0qbZt26Z8Pq81a9bosssu69PXYOvTpToaJskyyJ9LetXMfi3pTUn/5u6lWYQBgBS5+4hh2dzcrEcffVTvvvuuHnjggbPHQG/fvl0PPfSQOjs7VV9fr+7ubt1222164YUXdN5552nVqlXavXt3n75yuZzy+fzZ+/l8XhdffHFJnsuIYe3uB9y9rni70t23lWTPANBba2vhtmJF4dZzfxyampr09NNPq7u7W5IGXQY5ceKEZs+erTNnzmjXrl1nv//BBx9oyZIl2rp1q6qqqtTZ2akDBw7okksu0fr163X99ddr7969ffqaPXu2ZsyYoddff13urp07d+qGG24Y13PowenmACatK6+8Ups2bdKKFSs0bdo0LVq0SDt27OjT5sEHH9SSJUtUXV2tBQsW6MSJE5KkjRs3qqOjQ+6upqYm1dXV6eGHH9aTTz6piooKXXTRRbr//vsH7POxxx5Tc3OzPvnkE61evVqrV68uyXOxUh0D2FtDQ4Pz4QOTh1nhaIyslUsdSG7//v264oorsi6j7Aw2LmbWPtzZ4ZwUAwABENYAEABhDQABENYAEABhDQABENYAEABhDWDKSfMSqZs2bdLcuXN1wQUXlLRfwhoASui6667Tm2++WfJ+CWsAZaVxR6MadzSWrL+JvESqJF199dWaPXt2yervwenmACatib5EapoIawBloWc2vefQnj73W5tbx9znRF8iNU0sgwCYtCb6EqlpYmYNoCz0zKBLMaPu0dTUpBtvvFH33HOPZs6cqaNHjw6YXfe/ROqcOXMknbtE6pIlS/Tiiy+qs7NTx48fP3uJ1AMHDmjv3r1auXLluOtMgpk1gEmr9yVS6+rqdO+99w5o03OJ1GuuuUaXX3752e9v3LhRCxYsUG1trZYvX666ujo99dRTqq2tVX19vd5//32tXbt2QH/33XefcrmcTp06pVwupy1btpTkuXCJVIyoXC5NWi51IDkukTo4LpEKAJMUYQ0AARDWAFKVxlJrZGMdD8IaQGoqKyvV3d1NYBe5u7q7u1VZWTnqbTl0D0Bqcrmc8vm8urq6si6lbFRWViqXy416O8IaQGoqKio0f/78rMuYFFgGAYAACGsACICwBoAACGsACICwBoAACGsACICwBoAACGsACICwBoAAEoe1mU0zs7fN7KU0CwIADDSamfUGSfvTKgQAMLREYW1mOUlfl/TDdMsBAAwm6cz6B5Luk/THoRqY2TozazOztslwha2amsLHSGV5q6nJehQwmHJ4bfD6mHpGvOqemX1D0ofu3m5mjUO1c/cWSS1S4TMYS1ZhRg4dyv7z/syy3T8GVw6vDYnXx1STZGa9TNL1ZnZQ0k8lrTSzJ1OtCgDQx4hh7e7fc/ecu9dIulXSbnf/VuqVAQDO4jhrAAhgVJ8U4+6tklpTqQQAMCRm1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAEQ1gAQAGENAAGMGNZmVmlmb5rZr81sn5n940QUBgA4Z3qCNv8naaW7nzSzCkmvmtl/uPvrKdcGACgaMazd3SWdLN6tKN48zaIAAH0lWrM2s2lm9o6kDyX9wt3fSLcsAEBvicLa3T9393pJOUmLzay2fxszW2dmbWbW1tXVVeo6AfRTXS2ZZX+rqcl6JKaGUR0N4u7HJLVKunaQx1rcvcHdG2bNmlWi8gAM5eBByT3726FDWY/E1JDkaJBZZval4tfnSfqapPfTLgwAcE6So0FmS/oXM5umQrg/7e4vpVsWAKC3JEeD7JW0aAJqAQAMgTMYASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASAAwhoAAiCsASCAEcPazOaa2S/NbL+Z7TOzDRNRGADgnOkJ2nwm6bvu/paZzZDUbma/cPffpFwbAKBoxJm1ux9x97eKX5+QtF/SnLQLAwCck2RmfZaZ1UhaJOmNQR5bJ2mdJM2bN68EpaG6WjLLuopCHQCylTiszewCST+T9B13/7j/4+7eIqlFkhoaGrxkFU5hBw9mXQGAcpHoaBAzq1AhqHe5+3PplgQA6C/J0SAm6QlJ+939++mXBADoL8nMepmkOyStNLN3ire/TrkuAEAvI65Zu/urksrgz1wAMHVxBiMABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAhDUABEBYA0AAI4a1mf3IzD40s/cmoiAAwEBJZtY7JF2bch0AgGGMGNbu/itJRyegFgDAEKaXqiMzWydpnSTNmzevVN0CZ1VXS2ZZV1GoA+eUw/9LdbV08GC2NaStZGHt7i2SWiSpoaHBS9Uv0GOyvxmjKof/l6x/WEwEjgYBgAAIawAIIMmhez+R9J+SvmxmeTP72/TLAgD0NuKatbt/cyIKAQAMjWUQAAiAsAaAAAhrAAiAsAaAAAhrAAiAsAaAAAhrAAiAsAaAAAhrAAiAsAaAAAhrAAiAsAaAAEr24QOTSmOjpNZxbNtLa8J+erZL2r5U+x1N/++8I9XXD3yslPsa7ziMdnc7CvtrbR7b/say/VDbJO1rNPscbX1J2o93zMa7/6mKmTUABGDupf8EroaGBm9rayt5v6nrmdXt2SOTy1cU7yeZ5fXato8VK4bvo/92I7Uv1X6T6plRHz8+dJtS7Gu84zDa3RVncHsOFfa3orqwv9HOQEez/VDb9Bipr9Hsc7T1JWk/3jEbznj7NpNSiLIJZWbt7t4w1OPMrAEgAGbWg2lslO1pHdtPatasx7+fUvc53O5Ysx51+3Jcs2ZmDQAoC8yshzAZflIDU8VkeL8yswaASYCwBoAACGsACICwBoAACGsACICwBoAACGsACICwBoAACGsACICwBoAACGsACICwBoAACGsACICwBoAAEoW1mV1rZr81s9+Z2d+nXRQAoK8Rw9rMpkn6Z0mrJX1F0jfN7CtpFwYAOCfJzHqxpN+5+wF3/1TSTyXdkG5ZAIDepidoM0dSZ6/7eUlL+jcys3WS1hXvnjSz346/vGyZjXqTKkkflb6SkBiLvhiPvko+HmN4v5aLnrGoHq5RkrAebAgGfICOu7dIaklU2iRlZm3DfSzPVMJY9MV49MV4nJN0LJIsg+Qlze11PyfpD2MtDAAweknC+r8kXWZm883sTyTdKumFdMsCAPQ24jKIu39mZn8n6RVJ0yT9yN33pV5ZTFN6GagfxqIvxqMvxuOcRGNhHv3z2wFgCuAMRgAIgLAGgAAI6xLgdPxzzGyumf3SzPab2T4z25B1TVkzs2lm9raZvZR1LVkzsy+Z2bNm9n7xNbI065qyZGb3FN8n75nZT8yscqi2hPU4cTr+AJ9J+q67XyHpaknfnuLjIUkbJO3Puogy8U+SXnb3yyXVaQqPi5nNkbReUoO716pwAMetQ7UnrMeP0/F7cfcj7v5W8esTKrwZ52RbVXbMLCfp65J+mHUtWTOzP5W0XNITkuTun7r7sWyrytx0SeeZ2XRJX9Qw57AQ1uM32On4UzacejOzGkmLJL2RbSWZ+oGk+yT9MetCysAlkrok/bi4LPRDMzs/66Ky4u6HJT0i6feSjkg67u4/H6o9YT1+iU7Hn2rM7AJJP5P0HXf/OOt6smBm35D0obu3Z11LmZgu6auSHnP3RZL+V9KU/RuPmf2ZCr+Fz5d0saTzzexbQ7UnrMeP0/H7MbMKFYJ6l7s/l3U9GVom6XozO6jC8thKM3sy25IylZeUd/ee37SeVSG8p6qvSfpvd+9y9zOSnpP0F0M1JqzHj9PxezEzU2FNcr+7fz/rerLk7t9z95y716jwutjt7kPOnCY7d/8fSZ1m9uXit5ok/SbDkrL2e0lXm9kXi++bJg3zB9ckV93DMDgdf4Blku6Q9K6ZvVP83j+4+79nWBPKx92SdhUnNgck/U3G9WTG3d8ws2clvaXCUVRva5hTzzndHAACYBkEAAIgrAEgAMIaAAIgrAEgAMIaAAIgrAEgAMIaAAL4f5LMfsxVPfI8AAAAAElFTkSuQmCC)
%% Cell type:markdown id: tags:
### Exercize 2 : sklearn implementation
### Exercise 2 : sklearn implementation
- Compare the results obtained with the simple code above
- Comment and explain the role of the input parameters used in this implementation
%% Cell type:code id: tags:
``` python
#https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
# https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters = 2, init = 'k-means++', max_iter = 10, n_init = 10, random_state = 0)
kmeans = KMeans(n_clusters=2, init="k-means++", max_iter=10, n_init=10, random_state=0)
kmeans.fit(X)
y_kmeans = kmeans.fit_predict(X)
Y0=X[y_kmeans==0]
Y1=X[y_kmeans==1]
plt.scatter(Y0,np.zeros_like(Y0)+.7,c='red',marker='o', label='class 0 skl')
plt.scatter(Y1,np.zeros_like(Y1)+.7,c='green',marker='o',label='class 1 skl')
plt.scatter(X0,np.zeros_like(X0)+.5,c='red',marker='+', label='class 0')
plt.scatter(X1,np.zeros_like(X1)+.5,c='green',marker='+',label='class 1')
Y0 = X[y_kmeans == 0]
Y1 = X[y_kmeans == 1]
plt.scatter(Y0, np.zeros_like(Y0) + 0.7, c="red", marker="o", label="class 0 skl")
plt.scatter(Y1, np.zeros_like(Y1) + 0.7, c="green", marker="o", label="class 1 skl")
plt.scatter(X0, np.zeros_like(X0) + 0.5, c="red", marker="+", label="class 0")
plt.scatter(X1, np.zeros_like(X1) + 0.5, c="green", marker="+", label="class 1")
plt.legend()
```
%%%% Output: execute_result
......
......@@ -11,130 +11,133 @@
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import scipy.stats as stats
%matplotlib inline
```
%% Cell type:markdown id: tags:
## Create data set :
%% Cell type:code id: tags:
``` python
D1=np.random.randn(80,)*.1 +1
P1=np.random.rand(80,)*2*np.pi
D2=np.random.randn(40,)*.2
P2=np.random.rand(40,)*2*np.pi
C1=np.zeros((80,2))
C1[:,0]=D1*np.cos(P1)
C1[:,1]=D1*np.sin(P1)
C2=np.zeros((40,2))
C2[:,0]=D2*np.cos(P2)
C2[:,1]=D2*np.sin(P2)
D1 = np.random.randn(80,) * 0.1 + 1
P1 = np.random.rand(80,) * 2 * np.pi
D2 = np.random.randn(40,) * 0.2
P2 = np.random.rand(40,) * 2 * np.pi
C1 = np.zeros((80, 2))
C1[:, 0] = D1 * np.cos(P1)
C1[:, 1] = D1 * np.sin(P1)
C2 = np.zeros((40, 2))
C2[:, 0] = D2 * np.cos(P2)
C2[:, 1] = D2 * np.sin(P2)
plt.subplot(121)
fig=plt.scatter(C1[:,0],C1[:,1],marker='+', color='blue')
fig=plt.scatter(C2[:,0],C2[:,1],marker='o', color='red')
plt.axis('equal')
plt.title('theoretical')
X=np.append(C1,C2,axis=0)
fig = plt.scatter(C1[:, 0], C1[:, 1], marker="+", color="blue")
fig = plt.scatter(C2[:, 0], C2[:, 1], marker="o", color="red")
plt.axis("equal")
plt.title("theoretical")
X = np.append(C1, C2, axis=0)
plt.subplot(122)
plt.scatter(X[:,0],X[:,1])
plt.axis('equal')
plt.title('observed');
plt.scatter(X[:, 0], X[:, 1])
plt.axis("equal")
plt.title("observed");
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAArwUlEQVR4nO2df5QdVZXvP7ubBjqgaQJRoSEJo4jK8AjSD38wbwiIRmEQxF8ojiydWVk641rOLM2YGd5TnDcMmce88cegYnTx1CGK+IMYHtH4A4Ij83BMTKJEiEYkkA5KCOkopJFOst8fVZVUblfVrXvr16mq/VnrrntvVd2qU9Wnv+ecvffZR1QVwzAMo/kMVF0AwzAMoxxM8A3DMFqCCb5hGEZLMME3DMNoCSb4hmEYLcEE3zAMoyWY4BeIiMwTERWRw6ouSxQi8oSI/EHGc3xORP4hrzIZ9aeOdUJErhaRm6ouR9GY4OeMiDwoIhdUXY5ORGSNiPx5eJuqHq2qD1RVJsMwysUEvwG4OoIwjKKwOt8fJvg5IiL/BswBbhORJ4A3+buuEJGHROQxEbkqdPyAiCwRkV+KyE4RuUVEZoX2v1ZENonIhN9Df2Fo34Mi8gER+QnwpIgcJiIvFZH/8I/fKCIL/GOvAf4bcL1vxrne364i8jz/87CI/G8R2Soiu0XkByIy7O/7ioj82t/+fRE5rcDHaNQEEXmhXy8n/Hr62tDu40TkOyLyOxG5S0Tm+r8REfmIiDzq16efiMgf+vuOEJF/9v9XfiMiN4Tq4AIR2ebX+V8D/0dE7hORPwmV5zD/f+zF/vfI/wd/38l+uX4nIt8Bjiv8gbmAqtorxxfwIHCB/3keoMBngGHgDOD3wAv9/X8F3AOcCBwBfBr4kr/v+cCTwCuBIeBvgC3A4aHrbABO8s89CuwELsRryF/pf5/tH78G+POOsirwPP/zJ/xjRoFB4OXAEf6+dwLP8Mv4UWBD6ByfA/6h6udur9Lr+ZBfH/8OOBw4H/gdcKpfJ34H/LFfZz4G/MD/3UJgHTACCPBC4Hh/30eBlcAsv77dBlzr71sA7AX+yT/nMPBBYHmoTBcB9/ufu/0//D/gX/xz/bFf3puqfq6F/92qLkDTXjGCf2Jo/38Cl/uf7wNeEdp3PDAFHAb8D+CW0L4BYBxYELrOO0P7PwD8W0dZVgNX+p9jBd8/9yRwRor7G/F/N9P/boLfwhfeiPHXwEBo25eAq/06cXNo+9HAPrzOyfnAz4GXdvxW8Do4zw1texnwK//zAuBp4MjQ/uf5Qj3D/74c+KD/Ofb/AW8Uvhc4KrTvi20QfLODlcOvQ5/34P0DAMwFbhWR/aH9+4BnAycAW4ONqrpfRB7G67kEPBz6PBd4o4hcHNo2BNyZonzHAUcCv+zcISKDwDXAG4HZwP7Qb3anOLfRTE4AHlbVcN3dysH6eaBuquoTIvI4cIKq3uGbFD8BzBGRW4H349W/GcA6EQl+KnijzYAdqvpU6LxbROQ+4GIRuQ14LXCmvzvp/+EEYJeqPtlR9pN6fQh1w2z4+dNL+tGHgdeo6kjodaSqjgPb8Sot4Nk+8SrkeMy1Hsbr0YTPdZSqLk1RrseAp4DnRux7K3AJcAEwE2/UAt4/o9FetgMniUhYQ+ZwsH4eEE8RORrPTLMdQFU/rqpnAafhmS4X49XBSeC0UP2dqapHh84fVYe/BLwFr47+TFW3+NuT/h8eAY4RkaM6yt54TPDz5zdA2tj2G4BrQg6t2SJyib/vFuAiEXmFiAwB78Oz//9HzLluwuvpLBSRQRE50nd0nditXH4v7UbgX0TkBP/3LxORI/Bsqb/Hs3/OAP4x5b0ZzeaHeCaYvxGRId8hejFws7//QhH5IxE5HPifwA9V9WER+a8i8hK/Tj+J19HY59fBzwAfEZFnAYjIqIgs7FKOm4FXAe/GM8sExP4/qOpWYC3wYRE5XET+yC974zHBz59rgf8uIhPAG7oc+zE8J9W3ReR3eA7clwCo6mbgbcC/4vV+LgYuVtWno06kqg/j9XL+DtiB18NZzMG/8ceAN4jILhH5eMQp3g/8FPgR8Diec2wA+ALecHcc+JlfRqPl+PXwtcBr8OrnJ4G3q+r9/iFfBD6EV5fOAq7wtz8TT9h34dWrncA/+/s+gOcIvkdEfgt8F88JnFSOR/AcsC8Hvhza3u3/4a14/2uP++X8Qo+PoJaI77AwDMMwGo718A3DMFqCCb5hGEZLMME3DMNoCSb4hmEYLcHpiVfHHXeczps3r+piGA1l3bp1j6nq7LKva/XaKJKkeu204M+bN4+1a9dWXQyjoYjI1u5H5Y/Va6NIkuq1mXQMwzBaggm+YRhGS8hF8EXkRj+/9b0x+xf4ua83+K8P5nFdwygSq9dG08jLhv854HqSpyf/u6r+ScJ+w3CNz2H12mgQufTwVfX7eDkpDKMxWL02mkaZNvyX+cuMfTNpiTwRWSQia0Vk7Y4dO0osnmH0hdVrozaUJfg/Buaq6hl42R9XxB2oqstUdUxVx2bPLj1E2jB6weq1UStKEXxV/a2qPuF/XgUMiUg7Fg02GovVa6NulCL4IvIcf8UmRORs/7o7y7i2YRSF1WujbuQSpSMiX8JbZPg4EdmGt6DAEICq3oC3EMi7RWQv3jJml6sl4jccx+q10TRyEXxVfUuX/dfjhbcZRm2wem00DZtpaxiG0RJM8A3DMFqCCb5hGEZLMME3DMNoCSb4hmEYLcEE3zAMoyWY4BuGYbQEE3zDMIyWYIJvGIbREkzwDcMwWoIJvmEYRkswwTcMw2gJJviGYRgtwQTfMAyjJZjgG4ZhtAQTfMMwjJZggm8YhtESclnxynCfBQu89zVrqiyFYeTHivXjXLd6M9snJjlhZJjFC0/l0jNHqy6W05jgG6mwBsNwiRXrx/nbr/+Uyal9AIxPTPK3X/8pQE+i37ZGwwS/4YyMeO+7d3vvJtxGE7h65aYDYh8wObWP61ZvTi3YeTUadcIEPyeaKqTBfd1116Hfm3afxnRc7f2uWD/OxORU5L7tE5Opz3Pd6s2ZG426YYLfUAJhDnr2M2d67ybURhpc7v1et3pz7L4TRoZTnyeuceil0agbJvgZaXoPOLiPpt2XkYzLvd8kQV688NTU5zlhZJjxiHP10mjUDQvLbChr1nivc8/1XhMT3ssw0uBy7zdOkI+ZMdRTY7R44akMDw0esm14aLCnRqNuWA8/I3XoAW/Y4JUvS9lcvC+jOFzu/S5eeOoh5iYAAS76L8f3dJ6gcXDRT1EUufTwReRGEXlURO6N2S8i8nER2SIiPxGRF+dxXaM7a9bA/Pne5wULDjZMRnfaXK9d7v1eeuYorz9rFAltU+Cmex7izL//NivWj/d0rruXnM+vll7E3UvOb7TYQ349/M8B1wNfiNn/GuAU//US4FP+e2OI6gFX3evv9C8EjlsjNZ+jpfU6qffrQvTOnffvQCO279oz5Yxz2UVyEXxV/b6IzEs45BLgC6qqwD0iMiIix6vqI3lcv1eqFuKqsFj83qhbvc5biC89c3Ta712J3knyJbjiXHaRsmz4o8DDoe/b/G3T/jFEZBGwCGDOnDmlFC5vXIncCa7XOfnKyA1n6nVZQuxK9E6cjyHABeeyi5Ql+BKxLWpEhqouA5YBjI2NRR7TL64IcdkENvyApt9viThRr6E8IXYleifKcRvGBeeyi5Ql+NuAk0LfTwS2l3Tt0nEtcqezPFlx5b4cwJl6XZYQuxK9EzRiV6/cNG3Wbb/OZRd8E0VTluCvBN4jIjfjObV2V2HnTCPETRaztPfU2TA08VnkhBP1GpKFOE8hi+pZVxW9E/gY0txfcMz4xCSDIuxTZbTDEe2Cb6JochF8EfkSsAA4TkS2AR8ChgBU9QZgFXAhsAXYA7wjj+u6TpFCWUXD1DaTWJ3qdZwQn/eC2bkKWdGx6/00TlHO5c5zhp/BPvUsauFn4YpvomjyitJ5S5f9CvxlHtfKg6SefR5iVlch7HwGndvbRp3qdZwQ5ylknWL8kTfPz1UMi+plRz2DgOBZuOKbKBqbaVszquxlu+abMA4lqqf711/eEHlsr0JWhsmjqF52t3sNGjAXfBNFY4Lvk4eY1d3kEefcjSt/r/eX5vi6PTPXyUvIyjB5FNXL7hbCGYyGotI1nPeC2V3P380M5ZIz2AS/ZrjQy+4W9WOi7Q55OVnLMHkU1ctOCuEMnsWlZ46yduvjLL/noQNxtQp8bd04Y3NnAfGzjpNGPq45g03wO8gjwVjdBa9buXsdyaQ5vu6jI1fJy8lahsmjqAig8DOIi9KB6HQNk1P7+PBtm3hqan+kaHcb+bjmDDbBrykuCmGcaBvV0i2KJQ1lhGMWGQGU5hnEjVZ27Zm+ulZaZ69rzmAT/AJwUYzzJG4kEyfwaUY+TRkdNZWyUgnn0Tj1SzdbfydpnL2uOYNN8I3cMNFuNlWKcRnEjWKOOGwgcg3dOGdveOTj0kQ1MME3MtCroKc53hoJoyriRjFArGh3G/m4tsiKCb6RO/2Kto0MjKpJGsUkiXrcb1wKyQQT/NZgYmo0jTLFtB9zlmshmWCCH0ucQ9IEczpZn00vYZv2/LPRKZLnvWA2d96/w5keaFqKFtM8GpNuIZlV9P5N8BuOS/HtJtrVEiWSN93z0IH9LvRA01JkfHuvjUmccCeFZFbV+zfB76BTIDtXizLROkhejUnSKMqef34kJRELqEuGyCLj23tpTJKEOy4kc2TGEO+7ZeOBrJ3drpEnJvgNxwWTVNqGYcOGcsrTVtKKYR0yRBYZ395LY5LUOCxeeCqLv7qRqX0HhX1wQHjiqb3TxL7btfOi9YLfKT7dJhVZz/IgeT+b8FKMwTmDnv3Mmflcoy1EmRnSTiyqQ4bIIuPbe2lMujYOHbq+b7+SNMYq+tm3XvDbQllCGSX+3RqGqBHAhg3T1+I10hFnZnj9WaN8bd14olmnyklBvVBkfHtUYzI0KDz5+72cvOT2Q66V1Dhct3ozU/vTL19cxrNvreB3MzN0ipL1LOMp4tnMn++d10ZWvRNnZrjz/h1ce9npjYjSgeJm/nY2JiMzhnjiqb0HZtuG7fRJI424tQiiGBTh2stOtygdox6ksdPHibaZzfIlyczQ9PQIeRF+TucsvWNaArXATn/3kvOB6JFGkJ2zE+FQS8/w0GApYg8tFnwTmXpgf5fecS1hV93pZqePa0Tjev+vP2u0slFVawXfyJc8GlAT93xwLWFX3em3AXUtjw7AQGVXdoQ1a0xoXGHBAi/u3vLoZ+PSM0e59rLTGR0ZRoDRkeHSTAZNZPHCUxkeGjxkWy8N6JO/34vi2f4/fNsmVqwfL6CU6WhcD9/WTa0We6b5k3YKvmuJuppCvz31FevHWfyVjYdE6uzaM8Xir2485Lxl0jjBN+pHEIYZxNzfdZfX0w8iddpM2in4LibqahL9OLvjwjKn9qktcZgVWze1+bTx75V2mr9ra6caybNmw/vKHJk1RvCN+hJ2+AYTrtok6kmknebv2tqpRvKSiYHDt+yRWS6CLyKvBj4GDAKfVdWlHfsXAN8AfuVv+rqq/n0e1w6wdVObS5Ujs6rrdtoIEQvFdI/FC0+dZsMHb9Zu4PAte2SWOUpHRAaBTwCvAV4EvEVEXhRx6L+r6nz/lavYG81gzRqYmHCnIXahbqeNEMkaSWLkz6VnjnLdG89gZHjowLZjZgxx3RvOOCDmZY/M8ujhnw1sUdUHAETkZuAS4Gc5nLtnbN3U5lHhyKzyup02QsTFmG+ju7O37JFZHoI/Cjwc+r4NeEnEcS8TkY3AduD9qrop6mQisghYBDBnzpwcimcYfZNb3c5Sr9NGiFjahPpR9iS5PARfIrZ1xiL9GJirqk+IyIXACuCUqJOp6jJgGcDY2Fj6VHMxmL2+OVTwN8ytbudZry3evjmUPTLLQ/C3ASeFvp+I19M5gKr+NvR5lYh8UkSOU9XHcri+YRSFc3Xb4u2bR5kjszwE/0fAKSJyMjAOXA68NXyAiDwH+I2qqoicjecs3pnDtWOxmPvmU8Lf1Lm6bfH2zafIEVxmwVfVvSLyHmA1Xujajaq6SUTe5e+/AXgD8G4R2QtMAperxqzxZSSzfDlcdRU89BDMmQPXXANXXFF1qXIlq5Dn1RC4WLebHG/fdFNVmvtLGsFBdtNPLnH4qroKWNWx7YbQ5+uB6/O4VloaGXO/fDksWgR79njft271vkPjRD+JMkdvrtXtpsbbN91Ulfb+4kZwV6/cxO/37s/8fFqfLbNWXHXVQbEP2LPH294AFizwXnfd5b2C72X9vg40Nd4+yVTVBNLeX9xIbWJyKpfn0/jUCo3o2Qc89FBv2xtKt4Xmm0xT4+2bbKqC9PeXdqH5bueNo/GC3yjmzPHMOFHbG0BWM1wjzXgRNDHevqmmqoC09xcXl3/k0MC0ZRajft8NM+nUiWuugRkzDt02Y4a3vYXY4jXNoammqoC09xe3eM2HLj4tl+djPXwXiYvECRyzDY/SySri1gi4SVKUSlNNVQG93F/SCC7r8xGXoyPHxsZ07dq1VRejXDojccDrxS9b1jhhrxoRWaeqY2Vft431ujNKBbweqi29mD9J9brRJp1aRmnkEYmzfDnMmwcDA9778uV5ltComBXrxzln6R2cvOR2zll6R6VrpKYlryicOt67S5hJxzWyRuJYrH6jqWu8eh5ROHW9d5doZA+/1vHYcRE3aSNx3vveRsfqt526xqvHRZP0EmXy4ds21fLeXaKRgl9rskTiLF8OO2PSuLQsVr+p1DVePWsUzor145FhieD+vbtEI006tY7HzhKJk9SLb0isftupa7x61iicpF686/fuEo0U/NoTDsHshaRefEtj9ZtG2Qtm5EmWCWNJvfg63LsrmEmnScyaFb39qKPMYdsQLj1zlNefNcqgeGuzDIrw+rOaN/O2k5mhdWHDDA8NNP7e88QEvwrKDps88shiz2+Uxor143xt3Tj7/Pkz+1T52rpxJ8ITiwyZlKi1x4AjO/wCRjKNFHyno3SCsMmtW0H1YNhkkuinbSAefzx6+86dFo/fEFyN0glCJscnJlEOhkx2E/20jcREjMN2154pi8fvAbPhl03SxKoos0svcfWzZsVH6YR/F5SjwekZmoqrUTr9rMSVNq5+xfpxBkQOjGo6yXuRkCbT6NQKTkbpDAx4PftORGD//unb582LzpAJMHfuQWfse98bL/Zhjj0WJictdQP1TK1wztI7IqN0RkeGuXvJ+VmL1jcnL7l92uru4K0C/6ulF0X+Ju5ewLufxQtPZe3Wx1l+z0OR5+7kmBlDPDW1v/XpG1qbWsFJep1YlRR5s3UrvPOd8I53pBN78I6ziVm1xdWskv1MrEoalYxPTLL4qxu5KaXYg2fecdHc5RKNFnwn0+f2OrGqW/z800/DVLR9sydsYlYtiEufW3UPtp+GqFv8/NS+fKwPVZu7XMJs+GXT68Sqa66Znj2zX2bMgOHh6NGATcyqDS4ugNLPxKqoOQVpETik5z88NMgRhw0wMZl9kZAm0yjBd9JmH0UvE6vCDUScLT+OoSE44gh44gnv+/AwvOlN8PnPT7fh28QsZ0nKI+8SvTZEwbEfvm1TbNqEKAR4+XNncc8Du9inemAuwtjcWbWdlFYWtRf82oh8FgLRf8c7pptvDjvMc/h2bj/22IPiHrBzp/f9yith1SqL0nGUsMDPHB7iyaf3HjBvNDFD5FNTEcEKwNCggMLU/oN9+UDsf/zQ7mlzEcbmzuLay06vReNYFY2w4W/Y4HDcfV5cdVW0rX7vXnjmMz2BF/Eid266CR57zBP1KAftqlXw4INeVNA113jnttz5TtAZzz4xOTXNlt0kR2RUOGfA1D7l6CMPY2R46IC/4iNvns+DOycTQ0DvXnI+v1p6EYsXnsp1qzdb7vwQtRX88OSq3bs90W8c4QlXSeacnTu916xZh/bWu+XW72cSmFEoSQIYpu6OyGDCVVxYZsCuPVNMTE4xMmPoQG89zVyEfieCNZ3aCn4n8+fDzJlw7rmORuf0SqcYp2HnTnjb2+AZz4Cjj47/XeCgzWN1LSNX0gp5nR2RYTFOy649U/zVlzdw2ge/FRumGX4mrs5IrppcBF9EXi0im0Vki4gsidgvIvJxf/9PROTFWa8ZiPq55x4U+fnzs57VIaLEOC1PPAFPPpm8f/ny7KtrtYCy63YaIa+7IzLtKCaKJ5+O/934xOQB042rM5KrJrPgi8gg8AngNcCLgLeIyIs6DnsNcIr/WgR8Kut1o2hEzz6gSNHdudMbPcRl17QQTaCauh0Vzz40IBwzY8ipuPssFCm6gekmLrtmnUdGeZBHlM7ZwBZVfQBARG4GLgF+FjrmEuAL6uVxuEdERkTkeFV9JOvFGyPwncyZ03sYZi/s2eOFac6YUdsQzRIitEqv21kXCqkDcYu45MXk1D6OHBpgeGiwdiGaRYfg5mHSGQUeDn3f5m/r9RgARGSRiKwVkbU7duzIoXg1JWpGblyO2H7ZudPLoTN37sEInxbm1Ekgt7rdS70OR5rcveT8Rok9RI9icq7ZTOyZcnJGchJlOJrz6OFH/a06/SppjvE2qi4DloGXZKqXgoyMeO8TE738qkKWL4+fcRs1IzfvHn/QgDz4YL7nLZigZ3/XXYd+L6Cnn1vd7rder1g/ztUrNx2YQXrMjCE+dPFpTgsXJPdUo0Yxeff4Zw4POTkjOYl+Mo72Sh6Cvw04KfT9RGB7H8e0izRpjztn5CZlzuwH1fi0zAZUXLdXrB9n8Vc2HjLxaNeeKRZ/dSPg7sSrNGmPO8U4TYgmwIDAoEDMXK0D5D0YLoMyHM15mHR+BJwiIieLyOHA5cDKjmNWAm/3IxpeCuzOw34fMDLivXbv9l7Bd6fpJyQyyswTh4g3ASsw18RRw4icqAitguz4ldbt61ZvPkTsA6b2qdPhhf2EREaZeaJQhV/840V89M3zGU1wwMYtmOIy/WQc7ZXMgq+qe4H3AKuB+4BbVHWTiLxLRN7lH7YKeADYAnwG+Ius1609/YREXnHFQZt7N2bN8o4PZtTG/aaGETllzaSuum4n9excDi/sp6cazgKaxIAIK9aPH/BzxB1fx2ic814we5p9MG9Hcy5x+Kq6SlWfr6rPVdVr/G03qOoN/mdV1b/095+uqv2vahLBxIT3mjnTewXfnabXvPgBgYh3E/2dO+EvQtrTa1rmGlBGGG6VdTtJtFwWtH57qt1EHLy8OYu/svGAI9PV9QF6JVirODyeE8h9gfrGzLStHVkF+MILuxsqP/Wpg2kSwqODmkbkOL1WcQEsXngqQwPT/8ZDg+K0oGUV4aiebpip/crVKzcB7q4P0CtRZjAF7rw/30jF2mfLDON8rz5Mr3nxwyxf7mW9TJNy4b3vjXcCN4gmZk0NRKtuUTpZ5hJE9XSjCOe9r1s0ThRx5q5g9nBecfmNEvza0Sn6gcO2myj3knYh7dKHNSAQ8yaKexxNELNeyJJ2oc7EhaYKHNieR2psM+lUSb/ZKmsYWVMUbTPz1IEsE4hcdkYXSdxktM6RTtYEcI0T/Fr9w/ebrbKXyJpjj+29XI7TqJxJKQhSCdclr3uWTJVpndHHzIjOlVNXonwRcWatLI1i4wTfKcL57KMWF0mTrz7q92nj8Q8/HD72sf7K3gvd7rNASozJrwQX87p3a4DS5quPOkfaePwPXXxahjtIT5mNbWdKjSJCThtjwy9xun060sykjUuXMGdOut93W+f2xhuLd9KmKafRN2VMt++FNLNo4+zRgVClOcd1qzcnzrwt497TlLNIohZ5zxpyaj38okhjrkkKzUz6fTgHz2BMb2ju3HIE15FFVJrWsw9wLa97GnNNt7DMuHNcvXIT5yy9g7/+8gYARmJSHHebnJUXVS+iUkTIae16+HE9d+ciONLMpE0KzfzTP43+fdCDDkR2X0REQ5kTqmwRlVyISzbWrbdcNmkaoG5hmXHnmJicOhBuOT4xydCgMDQgh6SXKHNSlQuNbd5RWrUT/NqQZK4JExcbH/f7wcHokMzBQS+FQi/x/HmQ9j6NWJJMB0UM67OQtgFKEqq02TGn9inHzBhixuGHVbI2gGuNbR7UxqSTNvzOmaF91pm0cb+P6tGDJ/b793tpF8q0nTcwZUPZdLPTuzSTNI9UBmkds+AlQatqbYCmpG0IYz38osgykzbp93GO2qp61Fnv0+hqOnBp8lUeK3JFnWPP03vZFZHhssredBNXHxNNMz2/IsbGxnTt2kNzUUXZ6J2x25dBZ1QMeD3qmuXFKYpe6oKIrFPVsSLLE0VnvY7LBT8own7VRghNNzrNWuD1puuYF6cIeln6MKle18akY/g0IAlaErWaOJcTcSaOfarOxN4XjWumq7zJEs+f51yM2pl0onr2zsTel0WDk6D1S1Rd2LAB5s93vz50mg4GRNjXMfKuMva+LFwyXeVJ1nj+OB/P+27pfeWz2gm+0Uxa23j7hMXu5CW3Rx7T1jwzdSfr5Lm4v/s+1Z4ngtVa8J2LvTcqI1wXNmzwPu/efTCiK3yM6zQxHLDNZI3nTwpj7XXkZzZ8wwmanhOnF5oYDthmsq5V2y2MtZeRX617+AFtFYY20q233oTorSaGA7aZbpPnukXgBJ/fd8vGab4d6G3k1wjBN5pD3cS5KJrqwGwjSQ14Wodu8DnrrGsTfKMW9OPUtcbDcIW4BrwXh24eIz8TfMMwjIro1aGbdeRngm/UAovIMppI2RFZFqVjGIZREWVHZFkP36gV1rM3mkTZEVmZBF9EZgFfBuYBDwJvUtVdEcc9CPwO2AfsLTthVdgM0KtJwEwI7aQOdTsI5xufmGTQT8cwmkIweknEZRRPmRFZWU06S4DvqeopwPf873Gcp6rzq8hOaBh94HTdDifUAg7EZ3dLrOXiouhGeWQ16VwCLPA/fx5YA3wg4zlzozOUb2TEm24f3hfXc297bhfD7bodFc4XkDTd3rVF0Y1yydrDf7aqPgLgvz8r5jgFvi0i60RkUdIJRWSRiKwVkbU7duzIWDzD6Jtc63be9brbdPpew/0sMVs76NrDF5HvAs+J2HVVD9c5R1W3i8izgO+IyP2q+v2oA1V1GbAMvIUierjGNKJC+UZGDt3Xy2+NfHHg2T5fRO6N2J573c6zXkP3dWFHZgz19DtLzJYPrvtHuvbwVfUCVf3DiNc3gN+IyPEA/vujMefY7r8/CtwKnJ3fLaRjwwZPYHbv9l5tXGjDmMbP61q3uyXUilvIzhKzFUcd/CNZTTorgSv9z1cC3+g8QESOEpFnBJ+BVwFRvarCWLPGWwij399a7z5f0i5IXzFO1+1ghag4dk9OXx82/LumrixVJUn+EVfI6rRdCtwiIn8GPAS8EUBETgA+q6oXAs8GbhWR4HpfVNVvZbxuz5iJxugR5+v2pWeOHgjL7CTJRGOJ2YqhDv6RTIKvqjuBV0Rs3w5c6H9+ADgjy3XyJFgcw6iOOjS+danbUal3BTjvBbOrK1RLqYN/pBWpFcImg/nz+zfv5F0Ww8hC4CDsNCMo8LV146XajrMs0t0U6uAfaU1qhcBp20tcvcs90CZgz7V/OvOod5IUW593JEnWRbqbQh0Wrmm04HdOnqrSnGMTuYw8SZp4FRBlOy5CnG0y10Fc9480WvA7mT/fE/3585Nn2AbHmDgbrpLGERjYjsM9+gE/506YrOJcB2el4dFowY9yDlZlP6+DozJP2nKfVdFt4lVgO+7s0UetiQrZxLkOzsq8cH1iVTda4bQNExdXH44ND/LtzJwJ555rsfiGe0Q5CMV/D8fWpzH9QDZxroOzMg/qMLGqG43u4Qe4JNYulaUIzFdRDmkdhGl67lnFuQ7Oyjxogq+iFYKfhraZXIz6k8ZBGGduGRRhv2pu4uy6szIPmuCrMME3csUaTreImpg1PDRo6RT6oAm+ChN8DhUnWwnLaAJh5+LM4SGOHBpgYs9UYo++7g7JoolrPOvkqzDBNwrBGsLq6IzMmZicYnhokI+8eX6sgNvkqe40wVfRasHvx8FYpFPSRg1GHvTjXCzSIdmkkUPdfRWtFnzDfawR7J1+nItFOSRt5BBNVY1gqwU/jYOxc18RTsmmhDLWtdxNo5tzMUpsinJINiGUsUm5h1o38cqoBzVZJMVJkiZCxU0eOu8FswuZPFX3UMYiJltVuVBKq3v4AUk9+7hed5692LqHMjZlhNIUkpyL5yy9I1Js7rx/B9dednruZoa6hzIWMUKpshE0wTecpO6NYNXEOReTxKYIh2TdQxmLEOcqG0ET/BjiBKdIAcpyziqF0cS5PsSJzYAIJy+5PXcHYpZQRheie4oQ5yobQRN8w2ms8ciXKLGBg1k0i3Ag9jNycCW6pwhxrjKe3wS/C529V9fs1C6Vq+pnYXSnU2yKyI+fB65E9xQlzlXF85vgG0bLCIvNyUtujzym6igal6J76j7ZKowJfkpctVO7Wi6jHrgaReNqueqOxeG3jH7i2S0Gvrm4unhJr+VasX6cc5bewclLbuecpXekipPv5zd1x3r4PeJqD9rVchlu42pCsF7K1Y+D1xWncNmY4LcE1xLFdcNMVOXhqo06bblcSxYXhwthpplMOiLyRhHZJCL7RWQs4bhXi8hmEdkiIkuyXNPojplgsmN12z3iTDAuJYuLw5X1cLP28O8FLgM+HXeAiAwCnwBeCWwDfiQiK1X1ZxmvbfRAP87dKhzCDoWZWt12iCQTTD8O3rKdwq6EmWbq4avqfaraLePP2cAWVX1AVZ8GbgYuyXLdOlFmb9sSjuWH1e1kynZ4JglmP47nsp3VroSZlmHDHwUeDn3fBrwk7mARWQQsApgzZ06xJWsh/fSUy+xd1yzMNHXdblK9rsLh2S0HEPTmeC7bWe1KmGlXwReR7wLPidh1lap+I8U1JGKbRmzzdqguA5YBjI2NxR7nOlWYJmomli7wfBG5N2J77nW7KfUaqjFPdBPMfhzPZTqrXUki11XwVfWCjNfYBpwU+n4isD3jOZ3CBDZ/SnqWP1fVWIdsChpbt5MiSqowT7gimP3iSvhrGSadHwGniMjJwDhwOfDWEq5bKVX2tq3hKY1G1u1uJpsqzBOuCGYWXAh/zST4IvI64F+B2cDtIrJBVReKyAnAZ1X1QlXdKyLvAVYDg8CNqropc8kdwKGIEiNn2ly3u5lsquptuyCYdSeT4KvqrcCtEdu3AxeGvq8CVmW5Vl2pi/hbY3Uoba7b3Uw2deptuzDZySVspm0GzElqNJE0Jps69Lbbmj4hCUue1nIsdt/oxNWEar1S5WLhrmI9/Bywnr3RJOpksknClclOLmGC33LMLGVEUQeTTTdcmezkEmbSMQyjkTTFNJUn1sM3AOvZG82jKaapPDHBNwyjsTTBNJUnZtIxDMNoCSb4hmEYLcEE3zAMoyWY4BuGYbQEE3zDMIyWYIJvGIbREkzwDcMwWoIJvmEYRkswwTcMw2gJJviGYRgtwQTfMAyjJZjgG4ZhtAQTfMMwjJZggm8YhtESTPANwzBaggm+YRhGSzDBNwzDaAmZBF9E3igim0Rkv4iMJRz3oIj8VEQ2iMjaLNc0jDKwum00kaxLHN4LXAZ8OsWx56nqYxmvZxhlYXXbaByZBF9V7wMQkXxKYxiOYHXbaCJl2fAV+LaIrBORRSVd0zDKwOq2URu69vBF5LvAcyJ2XaWq30h5nXNUdbuIPAv4jojcr6rfj7neImARwJw5c1Ke3jD64vkicm/E9tzrttVrwwW6Cr6qXpD1Iqq63X9/VERuBc4GIgVfVZcBywDGxsY067UNI4Gfq2qsQzYNaeu21WvDBQo36YjIUSLyjOAz8Co8h5hh1Bqr20bdyBqW+ToR2Qa8DLhdRFb7208QkVX+Yc8GfiAiG4H/BG5X1W9lua5hFI3VbaOJZI3SuRW4NWL7duBC//MDwBlZrmMYZWN122giNtPWMAyjJZjgG4ZhtARRdTdgQER2AFsTDjkOcGGGo5XjUOpSjrmqOruswgSkqNdQn2fYljJAfcoRW6+dFvxuiMjarGF1Vg4rh4u4UnYXyuFCGZpSDjPpGIZhtAQTfMMwjJZQd8FfVnUBfKwch2LlyI4rZXehHC6UARpQjlrb8A3DMIz01L2HbxiGYaTEBN8wDKMl1Erwe1h27tUisllEtojIkgLKMUtEviMiv/Dfj4k5Lvfl77rdm3h83N//ExF5cR7X7aMcC0Rkt3/vG0TkgwWV40YReTQmzXFpzyMrVretbndco5h6raq1eQEvBE4F1gBjMccMAr8E/gA4HNgIvCjncvwvYIn/eQnwTzHHPQgcl+N1u94bXp6XbwICvBT4YQF/hzTlWAD83xLqxB8DLwbujdlf+PPI6T6sblvdLrxe16qHr6r3qermLoedDWxR1QdU9WngZuCSnItyCfB5//PngUtzPn8cae7tEuAL6nEPMCIix1dQjlJQb7GRxxMOKeN5ZMbqttXtMEXV61oJfkpGgYdD37f52/Lk2ar6CID//qyY4/Je/i7NvZVx/2mv8TIR2Sgi3xSR03IuQ1rKeB5lYXXb6nZAX88iU3rkIpDsSypGrTrdc+xpUjl6OE3qpR3TFitiW+e95XL/OZTjx3g5PZ4QkQuBFcApOZcjDWU8j1RY3U4uVsQ2q9vx9PUsnBN8zb6k4jbgpND3E4HteZZDRH4jIser6iP+MOrRmHOkXtoxJWnuLZf7z1oOVf1t6PMqEfmkiBynqmUnnyrjeaTC6nYiVrd7o69n0USTzo+AU0TkZBE5HLgcWJnzNVYCV/qfrwSm9c6kmOXv0tzbSuDtvhf/pcDuYIieI13LISLPERHxP5+NV9d25lyONJTxPMrC6rbV7YD+nkWRnua8X8Dr8Fq23wO/AVb7208AVnV4sH+O522/qoByHAt8D/iF/z6rsxx4Xv6N/mtTXuWIujfgXcC7/M8CfMLf/1NiIj5KKMd7/PveCNwDvLygcnwJeASY8uvGn1XxPKxuW92uQ7221AqGYRgtoYkmHcMwDCMCE3zDMIyWYIJvGIbREkzwDcMwWoIJvmEYRkswwTcMw2gJJviGYRgt4f8D03NAMuM+s1oAAAAASUVORK5CYII=)
%% Cell type:markdown id: tags:
### Question 6
### Exercise 6
- Briefly explain why usual Kmeans algorithm will fail to detect the classes above
- Is the Kernel approach the only possibily for this kind of clustering problem?
- Is the Kernel approach the only possibility for this kind of clustering problem?
%% Cell type:markdown id: tags:
### Exercice 7
### Exercise 7
- Propose a change of representation space to allow successfull Kmeans clustering in a 1D space. Implement it (use Kmeans_basic.ipynb example)
- Propose a change of representation space to allow successful Kmeans clustering in a 1D space. Implement it (use Kmeans_basic.ipynb example)
- Explain the role of parameter 'gamma' , then change it in Kernel Kmeans code below and comment your findings
- Compare the initialization of this algorithm with the type of initialization uszed in the previous studies of Kmeans.
- Compare the initialization of this algorithm with the type of initialization used in the previous studies of Kmeans.
%% Cell type:code id: tags:
``` python
#Kernel computation
N=X.shape[0]
Ker=np.zeros((N,N))
# Kernel computation
N = X.shape[0]
Ker = np.zeros((N, N))
gamma = 5
gamma=5
for i in range(0, N):
for j in range(0, N):
d = np.sum((X[i, :] - X[j, :]) ** 2)
Ker[i, j] = np.exp(-gamma * d)
for i in range(0,N):
for j in range(0,N):
d=np.sum((X[i,:]-X[j,:])**2)
Ker[i,j]=np.exp(-gamma*d)
# Init
import numpy.matlib
converged = 0;
converged = 0
# Kernel K-means is sensitive to initial conditions (as is Kmeans). Try altering
# this initialisation to see the effect.
# this initialisation to see the effect.
K = 2;
Z = np.matlib.repmat(np.array([1,0]),N,1);
perm=np.random.permutation(N)[0:np.intc(N/2)]
Z[perm,:]=[0,1]
K = 2
Z = np.matlib.repmat(np.array([1, 0]), N, 1)
perm = np.random.permutation(N)[0 : np.intc(N / 2)]
Z[perm, :] = [0, 1]
di=np.zeros((N,K))
count=0
di = np.zeros((N, K))
count = 0
while converged == 0:
count+=1
Nk=np.sum(Z,axis=0)
converged=1
for k in range(0,K):
Vk=Z[:,k].reshape(N,1)
di[:,k]=np.diag(Ker)\
-(2/Nk[k])*np.sum(np.matlib.repmat(Vk.transpose(),N,1)*Ker,axis=1)\
+(float(Nk[k])**(-2))*np.sum( np.sum( \
(Vk@Vk.transpose())*Ker,axis=0), axis=0 )
oldZ=np.copy(Z)
Z=np.zeros((N,K))
for i in range (0,N):
if di[i,0]<di[i,1]:
Z[i,:]=[1,0]
if Z[i,0]!=oldZ[i,0] :
converged=0
count += 1
Nk = np.sum(Z, axis=0)
converged = 1
for k in range(0, K):
Vk = Z[:, k].reshape(N, 1)
di[:, k] = (
np.diag(Ker)
- (2 / Nk[k]) * np.sum(np.matlib.repmat(Vk.transpose(), N, 1) * Ker, axis=1)
+ (float(Nk[k]) ** (-2))
* np.sum(np.sum((Vk @ Vk.transpose()) * Ker, axis=0), axis=0)
)
oldZ = np.copy(Z)
Z = np.zeros((N, K))
for i in range(0, N):
if di[i, 0] < di[i, 1]:
Z[i, :] = [1, 0]
if Z[i, 0] != oldZ[i, 0]:
converged = 0
else:
Z[i,:]=[0,1]
if Z[i,1]!=oldZ[i,1] :
converged=0
#visu
IndC0=np.where(Z[:,0]==1)[0]
IndC1=np.where(Z[:,1]==1)[0]
plt.scatter(X[IndC0,0],X[IndC0,1],color='green',marker='o')
plt.scatter(X[IndC1,0],X[IndC1,1],color='cyan',marker='o')
plt.axis('equal');
Z[i, :] = [0, 1]
if Z[i, 1] != oldZ[i, 1]:
converged = 0
# visu
IndC0 = np.where(Z[:, 0] == 1)[0]
IndC1 = np.where(Z[:, 1] == 1)[0]
plt.scatter(X[IndC0, 0], X[IndC0, 1], color="green", marker="o")
plt.scatter(X[IndC1, 0], X[IndC1, 1], color="cyan", marker="o")
plt.axis("equal")
print('converged in {} iterations'.format(count))
print("converged in {} iterations".format(count))
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbWklEQVR4nO3db4hc13kG8OeRIoGUBCVaybFje3ftIkKTihQzqE4DJSV/sEWDEtOAwzYRbWBRXIP8KQgEKfmwH5pPVorjsGnVKvE2Jh/qWGnXOI6huIS6eBSsyI7jWjW78iITrxRYo8hYsvT2w71jjWbvvXNn7r9z7nl+IHZ3ZjRzdCW999z3vOe9NDOIiEj7bWh6ACIiUg8FfBGRQCjgi4gEQgFfRCQQCvgiIoF4T9MDyLJjxw6bnp5uehgiIt44ceLEOTPbmfSc0wF/enoa3W636WGIiHiD5HLac0rpiIgEQgFfRCQQCvgiIoFQwBcRCYQCvohIIBTwRSqwAGAa0X+w6fhnkaYp4IuUbAHALIBlABZ/nUU5QV8nEilCAV+kZAcBXBx47CKAwwXft8oTiYRBAV+kRAsAzqc8d6bgex9GNScSCYcCvjihLamKrOA7WfC9004YRU8kEg4FfGlcm1IVWcF3ruB7p50wip5IJBwK+NK4NqUq0oLvBICZgu89B2DrwGNbUfxEIuFQwJfGtSlVkRaUj5Tw3jMA5gFMAWD8dR7FTyQSDqe7ZUoYJhGlcZIe900v+B7EtcXbLSW/vwK8jEszfGlcG1MVb/V9fx7+rklIu5QS8EkeJfkGyRdSnifJ75A8TfJXJO8o43OlHdqWqmjTmoS0S1kz/H8BcFfG83cD2BX/mgXwcEmfKy0xA2AJwNX4q6/BHkhfe1iG3yWn4r9SAr6ZPQPgdxkv2QfgBxZ5FsAHSN5UxmdLGHyq089ae/C55FT8V1cO/2YAr/X9vBI/Jo5zIdD6VqeftCbRT+kdaUpdAZ8Jj1niC8lZkl2S3dXV1YqHJVlcCbS+5cT71yTS+FhyKv6rK+CvALi17+dbAJxNeqGZzZtZx8w6O3cm3nhdauJKoPWxTr+3JpEW9H0sORX/1RXwjwP4alytcyeANTN7vabPljG5Emh9binQxpJT8VdZZZk/AvDfAD5CcoXk10geIHkgfskigFcBnAbwfQD3lfG5Ui1XAq3PQbPqklMX1ljEH2VV6XzZzG4ys01mdouZ/ZOZfc/Mvhc/b2b2t2b2B2a228y6ZXyuVCtPoK0j4Phepz9qyWnvmBLRVngi+di6ssYi/qBZ4tqpEzqdjnW7Ojc0aQFRzv4Mopn9HK4FrF7A6c/zb4Vfwdg1Sce0Z/DYTiO5JcUUohOLhInkCTPrJD2n1gqSKWt26sqibpXqTpkkHdOewWPryhqL+EMBX8ZWZcBxITfdRMpk2LHrf96VNRbxhwK+jK2qgONKbrqJK5hhx67/eZ8Xs6UZCvgytqoCjiupoiZSJlm7dAeP7eAGr424dpxGPTmOckXlwtWXjEcBX8ZWVfWMK7npJlImSUEcSD+2M7h2krgSPzbqFdEoV1SuXH3JeFSlI86ZhhvVJ75UIU1j+PHKqrbK8/tH+Sxplqp0xCuu5KZ9qf8fdkU0bFY+yhWVK1dfMh4FfHGOS4HWhz79w1JPw9ZERkldqTLIbwr44iQfAq0rhl0RDZuVj3JF5crVl4xHAV/Ec8OuiIbNyke5onLp6ktGp0VbkZbzZfFZyqFFW5GAaVYuPQr4IgEoa01Em6789p6mByAifhhMDfXKOwFdLfhCM3yRCrRxJuxKywsZnwK+lCYpyLUx8A3TdPuBqo65Nl35TwFfSpEU5P4GwF8jvL4rTc6EqzzZ5N10FeJJ3hcK+FKKpCB3CcDlgcdCSAE0ORMe92STJ0jnveWlmqu5SwFfSjFKMGt7CqDJ9gPjnGzyBums8s7eCeOvoDy/yxTwpRSjBLO2911psv3AOCebUa4Kkso7+08Yadp+kveFAr6UIinIbQawaeCxEPquNLnRaZyTzShXBUmpn6z78Pa0/STvCwV8yS0rz5sU5I4C+GeEucOzqeZv45xsRlmMTUr9ZM3sgTBO8r5QwJdc8uR5k4Kcul7Wb/CYA9kLsnmvCtJSPxuRLqSTvA8U8CUXbbrxU94TdZ6rgrTUzxUknzAegU7yrlHAl1y06cZPeU/Uea7E0lI/vRNEiKk73yjgSy6605GfyjxRZ6V+lLrzgwK+5KI7HfmpzBO12iz7TwFfctF/dj/pRC391B5ZcutV3Yg/en9fhxGlcSZxLQUzKrVH9p9m+JJJjbD8V1Z+PW0B+OC4A5PaKeAHZNTgrUZY0i9tofc89G/CFwr4gRgneKv2XvoN68czSFeH7lHAD8Q4wVu199JvlH48ujp0kwJ+IMYJ3qq9l34zACZSnhv8N6GrQzcp4AdinOCtkj4ZdAT5/k3o6tBNpQR8kneRfJnkaZKHEp7/FMk1ks/Hv75ZxudKfuMEb9Xey6C8/yZ0deimwnX4JDcCeAjAZwGsAHiO5HEz+/XAS//LzP6i6OfJeMatx1btvQzK829iDtfX7AO6OnRBGTP8PQBOm9mrZnYJwKMA9pXwvlIy9TuRuujq0E1lBPybAbzW9/NK/NigT5A8SfIJkh9LezOSsyS7JLurq6slDE9EmqAJhnvKCPhMeMwGfv4lgCkz+ziAfwDwk7Q3M7N5M+uYWWfnzp0lDE9ERIByAv4KgFv7fr4FwNn+F5jZm2Z2If5+EcAmkjtK+GwpSJtjRMJRRsB/DsAukreR3AzgXgDH+19A8kaSjL/fE3/u+RI+WwrQ5hiRsBSu0jGzd0jeD+BJRLe3PGpmL5I8ED//PQB/CeDrJN8B8BaAe81sMO0jNcvaHKN8q0j70OW42+l0rNvtNj2M1tqA9YstQLQoc7XmsUgYFlBOq2ZJR/KEmXWSntNO24Bpc0y6hVMLmH5wGhu+tQHTD05j4ZQSXYOq6L6qNaVqKeAHTK0Tki2cWsDsT2exvLYMg2F5bRmzP51V0O9TRfdVrSlVTymdwOkSe73pB6exvLa87vGpbVNYemCp/gE5aBpRQB40hajmPsmwFOI47ynrZaV0dIvDwKl1wnpn1pJbfKU9HqJxu68mBfReClEN16qnlI7IgMltyasYaY+HqIruq1pTqp4CvsiAuU/PYeum60PT1k1bMffp0Fc3rqmi+6rWlKqngC8yYGb3DOY/P4+pbVMgiKltU5j//Dxmdiv51TNuc7Ss/jpquFY9LdqKiLSI6vBFYqqvl5CpSkeC0auvv3g5qgbv1dcDULpGgqAZfsBC29V4+OnD7wb7nouXL+Lw08Vvra0rB/GBZviB6u1q7IW/3q5GoL2LZFXV1+vKQXyhGX6ghm1zb6Oq6usPPnGwsisHkTIp4AcqxF2NVdTXL5xawPm3km/toJ254hoF/ECFuKuxivr6rFm8duaKa5TDD9Qcrs/hA2HsapzZPVNqXj1rFq+dueIazfADpV2N5UibxU9smdCCrThHAT9gWdvcJZ+9u/aC4HWPbd20FUfuPtLQiETSKeAHLrRa/DItnFrAsZPHYH1d3gli/8f3a3YvTlIOP2Ah1uKXKWkjl8Gw+MpiQyMSyaYZfsDaVItfxk7XUd9DN0oR3yjgB6wttfhl3IN2nPdIW7DdvmW72iyIkxTwA9aWWvwyeuSM8x5JC7abN27Gm2+/mXriUM8daZJy+AFrSy1+GamVUd5j4dQCDj5xcN0OW4LYtGETfn/599c93n/iUM8daZJm+AFrSy1+GT1y0l5rsOtm4r3UT1I7BYOtC/Y9Z9bOVNqtUyQPBfzAtaEWv4weOUnv0dOflkkK2nlMbpvUIq80TgFfvFdGj5z+90jSm4kPC84TWyZSTz5VdesUyUv3tBUZsOFbG67bTNVDEJPbJrG8tpz4+7Zu2or5z88DwLsnh8ltk5j79Bxmds+s65vf/3uUw5eyZN3TVou2IgPSgnoveA8GbSCa2R+5+8i7gTspgPceSzoZiNRBKZ2WGmyZcB/UQiGPhVMLuHDpwrrHe2mZpPTRI/c8gnPfOPdu4M4qvZzZPYOlB5Zw9e+uYumBJQV7qZVSOi002DIhyVb4WZFTpaSUC7B+9j7qeyhtI3XKSuloht9CSS0TBvnaQqFKaRU479v8vtzBWqWX4jIF/BbKW+SnYsDrpVXgLK8tv5ui2fHtHdjx7R2pO2VHLb3UzlupkxZtW2gSUefLPK8LXa+2/szaGWzgBlyxK+teQ/DdRdz+DVe9+vxfnPkFFl9ZzHyPpNLLwfSPdt5K1TTDb6E5RDn6LD62UCjbYMO0pEANILFEs+fi5Yt4uPtw5nukbQJT+kfqVkrAJ3kXyZdJniZ5KOF5kvxO/PyvSN5RxudKsqSWCV+H/y0UypaWs9/IjQCwrjHaOCa2TGD/x/fj8NOH16VttPNW6lY4pUNyI4CHAHwWwAqA50geN7Nf973sbgC74l9/AuDh+KtUZAYK6MOkBdardhUTWyYS++WM6u0rb+PYyWOJaZusen+RKpQxw98D4LSZvWpmlwA8CmDfwGv2AfiBRZ4F8AGSN5Xw2SJjy+pnX0awB4ALly6kpm3K6AEkMooyAv7NAF7r+3klfmzU1wAASM6S7JLsrq6uljA8kWRpAbcOZ9bOlNIDyCe6f3LzyqjSSUp0Dq5y5XlN9KDZPKIUMzqdjru7wsR7aa0OvvJvX6n8s3tXFzO7Z1ob4Pvp/sluKGOGvwLg1r6fbwFwdozXSIXuQ3R2Z/z1vmaH44ykVgej5tAntkzgkXsewcSWidy/Z++uvaMO1Wttun+yz8oI+M8B2EXyNpKbAdwL4PjAa44D+GpcrXMngDUze72Ez5Yc7kO0St4rGLwS/9zWoF90M1NWb/wkvZ24R+4+su73pVX6LL6yONKYfNeW+yf7rnDAN7N3ANwP4EkALwH4sZm9SPIAyQPxyxYBvArgNIDvo72xxknzIz7uszJuaD6YW++VaabpVfsk5eTTavhDK71sy/2TfafmaQHIqiZ3929/PNMPTieWOk5tm8LSA0tjvWdaf/w8713FeHyU1NBPDfyqoeZpgUubn2bPW/1UxWamYTn9C5cupF5BqPTymi19309Awb4JCvgBmB3xcZ9VcRvBYTn982+dT00bhVZ6maQ3u+/f2fBWQ2MJnVI6gbgP0YzqCqKZ/SyA7zY6ompU1Y8+T5O1iS0TOPeNc2N/RltNI7mZ3xSApVpHEgaldATfBfAOopz9O2hnsAeqmVH3B/vJbZOpTdbOv3Ve7Y0TqELHHZrhi2RIumIgmLqIq1n+etNIn+HPIarFP4OoYmcOyusXpRm+yJiSOmpmVeyU1YOnTZLadW8FsBdRanEZ0ZVnb/etrpGqo4AvrVTWnaRCq5evQlK77nlEm3O0+7ZeCvjSOmVsvuoZp82CrDeDaIH2avx1BsrtN0EBPyChdCss805So7RZ2LxxM47cfWTkzwiVdt/WTwE/EL1aaF/zpaOkaMrcfJW3zcJGbsTRfUeDqq8vKi23H96WtPoo4AfC526Fo6ZoRtl8ledE0t9R86pdTXzvq3bVqWDvw9VcWm7fnaPYPgr4gfA5XzpqiiZvO4Nxcv1V7OQtm09Xc0m5famOAn4gfM6Xjpqiybv5atiJJGn270NvHJ+v5qRa2ngVCJ+7FVbVcTKtCyZB/PCeH6a2aADW3yXLpXTOBiR3QSWimbS0mzZeBWBYztbnfGlVs+qs9EzW7D/pLlku8flqTqqlgN8CeXO2vuZLq+o4mXUiqaLNcl1U/SJplNJpgWmoG+G4Bhuj9dIzvt+4ZAHqUROqrJSOAn4LKGdbvqraLItUTTn8llPOtny6cUn9fNg74DvN8FsgrQJnP6IGVbqsF9f5XEXmGs3wWy6pAmc/gGPwY/ON+Kfs2bj2DtRDAb8lBitw1HpWqlLFTt602qdlKL1TJgX8lvK5lYK4rYrZeNZ6k65Oy6OA31JayJWqVDGZSNo70E9Xp+VQwG8pbb6RqlQxmehfh0qjq9PiFPBbyudWCuK2qiYTvXWotKCvq9PiFPBbzNdWCuK2vJOJcSt5dHVanfc0PQAR8c8MsicQg3X1vYXX3u8d9t6AWkNUQTN8ESld0UoeXZ1WQwFfREqnsmA3KeCLSOlUFuwmBXwRKZ0WXt2kgC8ipVNZsJtUpSMilRhWySP10wxfRCQQhQI+ye0knyL5Svz1gymvWyJ5iuTzJNXg3nFFW9/qRhYibio6wz8E4Gkz2wXg6fjnNH9uZn+c1phfmtUL0gTwFYzf+raK1rkiUo6iAX8fovtsIP76hYLvJw3oD9LA+vvjjrJhRjeyEHFX0YD/ITN7HQDirzekvM4A/IzkCZKzKa8BAJCcJdkl2V1dXS04PMkjKUgPyrthRhtuRNw1tEqH5M8B3Jjw1CiTtk+a2VmSNwB4iuRvzOyZpBea2TyiCi50Oh13b7jbInmC8fac7zWJa1cKg4+LSLOGzvDN7DNm9kcJvx4H8FuSNwFA/PWNlPc4G399A8BjAPaU90eQovIE4zeRLw+vDTftocX39ima0jmO6H7ZiL8+PvgCku8l+f7e9wA+B+CFgp8rJRp2tyEAuIx8l3TacNMOWnxvJ5qNnzUhOQHgx4gmiWcAfMnMfkfywwD+0cz2krwd0aweiFJI/2pmuSZ8nU7Hul1VcdZhAdfa0ab9iyCi7oXSftNITs1NIepeKe4ieSKtGrJQwK+aAn4zpqH/7KHbgOQTv0767ssK+NppK+soDy/qdtlOCviyjvLwopN+Oyngi8g6Oum3k7plCoDrF223IyrDvBw/N8r9SKU91O2yfTTDl3UleOdxLdj35G2PUHfttmrFRfLTDF9ytVYAhu/I7Z04eu9V9ZVB3Z8n4jvN8CV3n5thFRp1N05TozaR0SjgS65Su7QKjf6USlLtPlBd4zQ1ahMZjQK+JJbgbQYwgewKjcHcf5qqardVKz4arXeIAr4kluAdBXAO0a7KJSTnxPPk/qus3VateH7qjSOAAr7EZhAF9qwAPygrdVJH7bZqxfPTeocAqtKRAtJ639fZc0e14vlovUMAzfClAKVU/KH1DgEU8KUApVT8oZOzAAr4MqLBSg9g9Ny/1E8nZwEU8GUEZVV6qDywGeMszEu7KOBLbmVUeqg8UKQ5CviSWxmVHk2UB+qKQiSigC+5lVHpUXd5oK4o6qMTq/sU8CW3Mio96i4P1IajeujE6gcFfMktb6VH1kyv7vLA0Dcc1TXr1onVD9ppKyMZtrN1WI/63u/t3V1rElGwr6piJG03cAgbjuq8X0DoJ1ZfaIYvpcoz06uzPDDkDUd1zrq1k9cPCvhSKtdmeiFvOKrz7yLkE6tPFPClVC7O9ELdcFTn30XIJ1afKOBLqbJmem0o2/Ppz1D3rDvUE6tPFPClVGkzPcD/sj3fSg8165ZBNMu6OV2zOp2OdbvdpochJZhG873zi5qG/38GaT+SJ8ysk/ScZvhSC9cWc8fRhj+DhE0BX2rh4mLuqNrwZ5CwKeBLLdpQtteGP4OETQFfatGGBcQ2/BkkbAr4Ups2lO2N+2coWs7pUzmouEu9dEQqVrSnTZ09caTdNMMXqVjRnjYud6LUlYdfCgV8kl8i+SLJqyQT6z7j191F8mWSp0keKvKZIi7JE/CKlnO6Wg7q20Y0KT7DfwHAPQCeSXsByY0AHgJwN4CPAvgyyY8W/FyRxuUNeEXLOV0tB3X5ykOSFQr4ZvaSmb085GV7AJw2s1fN7BKARwHsK/K5Ii7IG/CKlnO6Wg7q6pWHpKsjh38zgNf6fl6JH0tEcpZkl2R3dXW18sGJH1zMFecNeEXLOV0tB3X1ykPSDa3SIflzADcmPHXYzB7P8RlMeCy1gY+ZzSPut9XpdNxt9CO1cbVKZZS7aQ27U9gwRX9/FeZw/d8L4MaVh6QbGvDN7DMFP2MFwK19P98C4GzB95SAZKVOmgyCoQe8um9XKcXVkdJ5DsAukreR3AzgXgDHa/hcaQlXc8Wuplrq1IbNdCEpWpb5RZIrAD4B4D9IPhk//mGSiwBgZu8AuB/AkwBeAvBjM3ux2LClTYbl513OFSvgiU8K7bQ1s8cAPJbw+FkAe/t+XgSwWOSzpJ3y5OdDT52IlEU7baVReUoblToRKYd66UijRiltVIAXKUYzfGmUy/n5Orm4z0DaRwFfGuXqLtI6qSeN1EUBXxql/Lx60kh9lMOXxoWen3d1n4G0j2b4Ig3TOobURQFfpGFax5C6KOCLNEzrGFIX5fBFHBD6OobUQzN8EZFAKOCLiARCAV9EJBAK+CIigVDAFxEJBM3cvW0syVUk3za0CTsAnGt6ECPQeKul8VbPtzG7Mt4pM9uZ9ITTAd8lJLtm1ml6HHlpvNXSeKvn25h9GK9SOiIigVDAFxEJhAJ+fvNND2BEGm+1NN7q+TZm58erHL6ISCA0wxcRCYQCvohIIBTwU5D8EskXSV4lmVpqRXKJ5CmSz5Ps1jnGgXHkHe9dJF8meZrkoTrHODCO7SSfIvlK/PWDKa9r9PgOO16MfCd+/lck76h7jAPjGTbeT5Fci4/n8yS/2cQ4+8ZzlOQbJF9Ied614ztsvE4d33XMTL8SfgH4QwAfAfCfADoZr1sCsMOH8QLYCOD/ANwOYDOAkwA+2tB4vw3gUPz9IQB/79rxzXO8AOwF8ASiVvZ3AvifBv8N5BnvpwD8e1NjTBjznwG4A8ALKc87c3xzjtep4zv4SzP8FGb2kpm93PQ48so53j0ATpvZq2Z2CcCjAPZVP7pE+wAci78/BuALDY0jS57jtQ/ADyzyLIAPkLyp7oHGXPr7zcXMngHwu4yXuHR884zXaQr4xRmAn5E8QXK26cEMcTOA1/p+Xokfa8KHzOx1AIi/3pDyuiaPb57j5dIxzTuWT5A8SfIJkh+rZ2hjc+n45uXs8Q36jlckfw7gxoSnDpvZ4znf5pNmdpbkDQCeIvmbeBZQuhLGy4THKqvLzRrvCG9T2/FNkOd41XpMh8gzll8i6rVygeReAD8BsKvqgRXg0vHNw+njG3TAN7PPlPAeZ+Ovb5B8DNFldSUBqYTxrgC4te/nWwCcLfieqbLGS/K3JG8ys9fjS/Q3Ut6jtuObIM/xqvWYDjF0LGb2Zt/3iyS/S3KHmbnQ9CuJS8d3KNePr1I6BZB8L8n3974H8DkAiav3jngOwC6St5HcDOBeAMcbGstxAPvj7/cDWHeF4sDxzXO8jgP4alxNcieAtV6qqgFDx0vyRpKMv9+DKAacr32k+bl0fIdy/vg2vWrs6i8AX0Q0u3gbwG8BPBk//mEAi/H3tyOqhDgJ4EVEqRVnxxv/vBfA/yKq5mhyvBMAngbwSvx1u4vHN+l4ATgA4ED8PQE8FD9/ChkVXY6M9/74WJ4E8CyAP214vD8C8DqAy/G/3685fnyHjdep4zv4S60VREQCoZSOiEggFPBFRAKhgC8iEggFfBGRQCjgi4gEQgFfRCQQCvgiIoH4fz0QAO4KaOsBAAAAAElFTkSuQmCC)
......
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment