#!/usr/bin/env python
# coding: utf-8

# <h1>Pau Machine Learning : PCA algorithm

# Voici le premier épisode de Pau ML, dont le but est d'échanger autour du data science. Pour commencer, nous attaquons une méthode très classique de data mining : l'analyse en composantes principales, ou PCA in english.

# In[1]:


get_ipython().run_line_magic('pylab', '--no-import-all inline')
from sklearn.decomposition import PCA
matplotlib.rcParams['figure.figsize'] = 10, 10


# <h2>1. Un exemple jouet

# In[2]:


alpha = 4
size = 100
a = np.random.uniform(-1, alpha - 1, (size, 2))
b = np.random.uniform(-2, alpha - 2, (size, 2))
c = np.random.uniform(-3, alpha - 3, (size, 2))
d = np.random.uniform(0, alpha, (size, 2))
e = np.random.uniform(1, 1 + alpha, (size, 2))
f = np.random.uniform(2, 2 + alpha, (size, 2))
X = np.concatenate((a,b,c,d,e,f),axis=0)
pca = PCA(n_components=2)
plt.scatter(X[:,0],X[:,1])
pca.fit(X)
print(pca.explained_variance_ratio_)


# La sortie ci-dessus renseigne le pourcentage de variance expliqué par chacun des axes de l'ACP. Ces valeurs sont calculées grâce à la diagonalisation de la matrice ${}^T\!XX$, qui représentent les corrélations entre les colonnes de $X$ (les variables de notre problème d'analyse de données). Chaque coefficient de cette matrice est défini par :
# $$
# ({}^T\!XX)_{ij} = \sum_{k=1}^n X_{ki}X_{kj}, \forall i,j=1, \ldots, p. 
# $$
# Mathématiquement, l'analyse en composantes principales est un simple changement de représentation (un changement de base): passer d'une représentation canonique (avec la base canonique) à une représentation idéale grâce à la matrice des corrélations (la base des composantes principales). Ci-dessous par exemple voici les nouvelles coordonnées de $X$ dans la nouvelle base :

# In[3]:


newX = pca.transform(X)
print(newX)
plt.scatter(newX[:,1], newX[:,1])


# <H2>2. Exemple digits

# L'intérêt de l'analyse en composantes principales est de considérer des données de dimension *plus grande que 2*. Comment représenter des observations ayant un grand nombre de feature ($p=10$, $p=100$ ou $p=10^6$) ? Comment réduire la dimension en tenant compte des corrélations entre les caractères ?  

# <H3>2.1 Description de l'échantillon

# Grâce au module sklearn, on va étudier un exemple concret : la reconnaissance d'écriture, et plus particulièrement la reconnaissance de chiffres. Dans cet exemple, chaque image représente un chiffre écrit à la main, à faible résolution (8*8 pixels). Nous sommes donc dans le cas où $p=64$.

# In[4]:


from sklearn import datasets
#iris = datasets.load_iris()
digits = datasets.load_digits()


# Voici par exemple une écriture du chiffre $0$ dans une matrice 8*8 pixels représentant les niveaux de gris de chaque pixel :

# In[5]:


digits.images[0]


# Voici une série d'images correspondant aux 10 premières observations du dataset digits :

# In[6]:


#Load the digits dataset
digits = datasets.load_digits()
#Display the first digit
for elt in range(10):
    plt.figure(1, figsize=(3, 3))
    plt.imshow(digits.images[elt], cmap=plt.cm.gray_r, interpolation='nearest')
    plt.show()


# <h3>2.2. PCA des digits 0 et 7

# Attaquons maintenant notre analyse en composantes principales en essayant de représenter notre échantillon avec seulement deux composantes principales (au lieu de 64 dimensions initiales).

# In[7]:


X = digits.data
pca = PCA(n_components=2)
pca.fit(X)
print(pca.explained_variance_ratio_)


# Le résultat ci-dessus nous dit que 27% de l'information est représenté par les deux premiers axes de l'ACP. Commençons par étudier le cas des digits 0 et 7. Comment sont représentés les observations de  0  et de 7 dans les deux premiers axes de l'ACP ?

# In[8]:


liste_0 = []
liste_7 = []
for k,elt in enumerate(digits.target):
    if (elt==0):
        liste_0.append(k)
    if (elt==7):
        liste_7.append(k)


# In[9]:


print('On étudie '+str(len(liste_0))+' écritures du chiffre 0 et '+str(len(liste_7))+' du chiffre 7.')


# In[10]:


X0 = digits.data[liste_0,:]
X7 = digits.data[liste_7,:]
X07 = np.concatenate((X0,X7))
pca = PCA(n_components=5)
pca.fit(X07)


# On va représenter ces écritures selon les deux premiers axes de l'analyse en composantes principales.

# In[11]:


Z07 = pca.transform(X07)
x = Z07[:,0]
y = Z07[:,1]
colors = list(digits.target[liste_0]) + list(digits.target[liste_7])
#area = np.pi * (15 * np.random.rand(N))**2
plt.scatter(x, y, c=colors, alpha=0.5)


# Not that bad ! Grâce à un simple changement de représentation, il est maintenant très facile de détecter automatiquement un 0 ou un 7 avec un algorithme de classification en deux dimensions comme une Analyse Discrimnante Linéaire (LDA in english). Si l'on regarde de plus près cette nouvelle représentation, voici à quoi ressemble le premier axe de cette représentation :

# In[12]:


components = pca.components_
components[0].reshape((8,8))


# Avec une image c'est mieux :

# In[13]:


plt.figure(1, figsize=(3, 3))
plt.imshow(components[1].reshape((8,8)), cmap=plt.cm.gray_r, interpolation='nearest')
plt.show()


#  Bien sûr, les choses se compliquent avec les 10 chiffres...

# In[14]:


X = digits.data
pca = PCA(n_components=5)
pca.fit(X)
Z = pca.transform(X)
x = Z[:,0]
y = Z[:,1]
colors = digits.target
#area = np.pi * (15 * np.random.rand(N))**2

plt.scatter(x, y, c=colors, alpha=0.5)
plt.show()


# # A vous de jouer !

# On étudie 84 peintures de Rembrant et Van Gogh grâce à l'ACP. Pour cela on représente chaque image par son histogramme de couleurs. Dans le cas présent, cela consiste à partitioner l’espace des couleurs (ici, $[0,255]^3$)
# en k parties égales et, pour chaque image I, calculer la proportion $p_{Ij}$ des pixels se trouvant dans la partie j, pour $j = 1, . . . , k$. Après avoir effectué cette étape, on associe à chaque image I le vecteur numérique de taille k contenant les proportions ($p_{I1},\ldots, p_{Ik}$).
# 
# ##Echauffement
# Les fichiers <a href="http://www.math.univ-angers.fr/~loustau/painting8.txt">painting8.txt</a> et <a href="http://www.math.univ-angers.fr/~loustau/painting64.txt">painting64.txt</a> contiennent l'histogramme des couleurs de chaque pour $k=8$ et $k=64$. Les 40 premières lignes correspondent aux tableaux de Rembrant et les 44 dernières à ceux de Van Gogh et les peintures se trouvent <a href = "http://www.math.univ-angers.fr/~loustau/painting.png"> ici </a>.
# 
# Faites une analyse en composantes principales pour séparer les peintures de Rembrant et Van Gogh.

# In[15]:


#PCA for painting dataset
paint_data = np.genfromtxt("painting64.txt", delimiter=',')
print(paint_data.shape)


# In[16]:


# PCA analyses
pca = PCA(n_components=3)
pca.fit(paint_data)
new_paint_data = pca.transform(paint_data)
x, y, z = new_paint_data[:,[0, 1, 2]].transpose()
print(pca.explained_variance_ratio_)

# plot
colors = 40*["blue"] + 44*["red"]
plt.subplot(211)
plt.scatter(x, y, color=colors)
plt.xlabel("First component")
plt.ylabel("Second component")
plt.subplot(212)
plt.scatter(x, z, color=colors)
plt.xlabel("First component")
plt.ylabel("Third component")


# Visualisation avec panda.

# In[20]:


import pandas as pd
from pandas.tools.plotting import scatter_matrix

df = pd.DataFrame({"Comp1": x, "Comp2": y, "Comp3": z})
pltmat = scatter_matrix(df, alpha=.6, figsize=(10, 10), diagonal='kde', marker="o")