<jupyter_start><jupyter_text>TME HMM : décodage des lettresLe but est le même que dans le TME précédent: apprendre à classer les lettres les lettres manuscrites enregistrées par un stylo intelligent.Les étapes sont un peu les mêmes que dans le TME précédent... Avec un modèle plus complexe.<jupyter_code>import numpy as np
import pickle as pkl
import matplotlib.pyplot as plt

with open('ressources/lettres.pkl', 'rb') as f:
    data = pkl.load(f, encoding='latin1') 
X = np.array(data.get('letters')) # récupération des données sur les lettres
Y = np.array(data.get('labels')) # récupération des étiquettes associées 

nCl = 26
# affichage d'une lettre (cf TME précédent)
def tracerLettre(let):
    a = -let*np.pi/180; # conversion en rad
    coord = np.array([[0, 0]]); # point initial
    for i in range(len(a)):
        x = np.array([[1, 0]]);
        rot = np.array([[np.cos(a[i]), -np.sin(a[i])],[ np.sin(a[i]),np.cos(a[i])]])
        xr = x.dot(rot) # application de la rotation
        coord = np.vstack((coord,xr+coord[-1,:]))
    plt.figure()
    plt.plot(coord[:,0],coord[:,1])
    #plt.savefig("exlettre.png")
    return
# discrétisation (cf TME précédent)
def discretise(x, d):
    intervalle = 360 / d
    xd = [np.floor(x[i] / intervalle) for i in range(len(x))]
    return np.array(xd)<jupyter_output><empty_output><jupyter_text>A. Apprentissage d'un modèle connaissant les états cachés du système A1. Hypothèse gauche-droiteOn fait l'hypothèse que les états sont connus... Alors que ce n'est pas le cas. Mais il existe des stratégies simples (et parfois efficaces) pour attribuer arbitrairement des états sur les chaines.La plus classique est l'hypothèse gauche-droite, qui est bien adaptée aux signaux courts et non périodiques:* On définit le nombre d'états N* On découpe chaque série d'observations en N portions à peu près égales* On affecte l'état 0 au début, puis on incrémente jusqu'à l'état N pour la dernière portion de la chaineSur un exemple:```X0 = [ 1.  9.  8.  8.  8.  8.  8.  9.  3.  4.  5.  6.  6.  6.  7.  7.  8.  9.  0.  0.  0.  1.  1.]S0 = [ 0.  0.  0.  0.  0.  0.  1.  1.  1.  1.  1.  1.  2.  2.  2.  2.  2.  3.  3.  3.  3.  3.  3.]```Au niveau de la mise en oeuvre, vous définirez la méthode ```def initGD(X,N):```qui prend un ensemble de séquences d'observations et qui retourne l'ensemble des séquences d'états. Pour chaque séquence $x$ vous pouvez utiliser:```np.floor(np.linspace(0,N-.00000001,len(x)))```<jupyter_code>def initGD(x,N):
    return np.floor(np.linspace(0,N-0.00000001,len(x)))
    
# 14 observations aléatoires
x = np.floor(np.random.rand(14)*10)
print(x)
print(initGD(x,4))
# Construire la base Q des états, initialisés en Gauche-Droite correspondant aux observations X
N = 5 # nombre d'états

Q = np.array([initGD(X[i], N) for i in range(len(X))])

print(X[0],"\n",Q[0])<jupyter_output>[ 36.214493 347.719116 322.088898 312.230957 314.851013 315.487213
 313.556702 326.534973 141.288971 167.606689 199.321594 217.911087
 226.443298 235.002472 252.354492 270.045654 291.665161 350.934723
  17.892815  20.281025  28.207161  43.883423  53.459026] 
 [0. 0. 0. 0. 0. 1. 1. 1. 1. 2. 2. 2. 2. 2. 3. 3. 3. 3. 4. 4. 4. 4. 4.]<jupyter_text>Check:```[ 36.214493 347.719116 322.088898 312.230957 314.851013 315.487213 313.556702 326.534973 141.288971 167.606689 199.321594 217.911087 226.443298 235.002472 252.354492 270.045654 291.665161 350.934723  17.892815  20.281025  28.207161  43.883423  53.459026] [0. 0. 0. 0. 0. 1. 1. 1. 1. 2. 2. 2. 2. 2. 3. 3. 3. 3. 4. 4. 4. 4. 4.]```  A2. ApprentissageEtant donné la structure d'un MMC:* les observations n'influencent pas les états: les matrices $\Pi, A$ s'obtiennent comme dans un modèle de Markov simple * chaque observation ne dépend que de l'état courantLa nature des données nous pousse à considérer des lois de probabilités discrètes quelconques pour les émissions. L'idée est donc de procéder par comptage en définissant la matrice $B$ comme suit:* $K$ colonnes (nombre d'observations), $N$ lignes (nombre d'états)* Chaque ligne correspond à une loi d'émission pour un état (ie, chaque ligne somme à 1)Ce qui donne l'algorithme:1. $b_{ij}$ = comptage des émissions depuis l'état $s_i$ vers l'observation $x_j$1. normalisation des lignes de $B$Donner le code de la fonction ```def learnHMM(allX, allS, N, K):``` qui apprend un modèle à partir d'un ensemble de couples (seq. d'observations, seq. d'états)***Variante stabilisée***: afin d'éviter les transitions à probabilité nulle et de régulariser l'ensemble du système, il est intéressant d'initialiser les matrices de transition à une valeur non nulle. Cette initialisation *spéciale* peut être faite de manière optionnelle en utilisant les arguments par défaut de python:<jupyter_code>def learnHMM(allx, allq, N, K, initTo0=False, eps=1e-5):
    if initTo0:
        A = np.zeros((N,N))
        B = np.zeros((N,K))
        Pi = np.zeros(N)
    else:
        A = np.zeros((N,N)) # np.ones((N,N))*eps => non pertinent pour un modèle GD
        B = np.ones((N,K))*eps
        Pi = np.ones(N)*eps 
                
    for i in range(len(allq)):
        for k in range(N):
            if(allq[i][0] == k):
                Pi[k] += 1
        
        n = len(allq[i])
        for j in range(n-1):
            A[int(allq[i][j])][int(allq[i][j+1])] += 1
            B[int(allq[i][j])][int(allx[i][j  ])] += 1
            
    B[int(allq[i][j])][int(allx[i][j])] += 1
                
    t = Pi.cumsum()
    Pi = Pi/t[N-1]
    A = A / np.maximum(A.sum(1).reshape(N,1),1)
    B = B / np.maximum(B.sum(1).reshape(N,1),1)
    
    return Pi , A, B
# le nombre d'états a été fixé au dessus (à 5)
K = 10 # discrétisation (=10 observations possibles)
Xd = discretise(X,K)
Pi, A, B = learnHMM(Xd[Y=='a'],Q[Y=='a'],N,K, True)

print(Pi,"\n", np.around(A,decimals=2), "\n", np.around(B,decimals=2))<jupyter_output>[1. 0. 0. 0. 0.] 
 [[0.79 0.21 0.   0.   0.  ]
 [0.   0.76 0.24 0.   0.  ]
 [0.   0.   0.77 0.23 0.  ]
 [0.   0.   0.   0.76 0.24]
 [0.   0.   0.   0.   1.  ]] 
 [[0.06 0.02 0.   0.   0.   0.   0.   0.04 0.49 0.4 ]
 [0.   0.04 0.   0.13 0.09 0.13 0.02 0.09 0.41 0.09]
 [0.   0.   0.   0.02 0.12 0.5  0.31 0.04 0.   0.  ]
 [0.07 0.   0.   0.   0.   0.   0.26 0.33 0.2  0.15]
 [0.71 0.1  0.   0.   0.   0.   0.   0.02 0.02 0.14]]<jupyter_text>Validation sur les séquences de la classe 'a':```  Pi=[ 1.  0.  0.  0.  0.]  A=[[0.79 0.21 0.   0.   0.  ]     [0.   0.76 0.24 0.   0.  ]     [0.   0.   0.77 0.23 0.  ]     [0.   0.   0.   0.76 0.24]     [0.   0.   0.   0.   1.  ]]   B=[[0.06 0.02 0.   0.   0.   0.   0.   0.04 0.49 0.4 ]     [0.   0.04 0.   0.13 0.09 0.13 0.02 0.09 0.41 0.09]     [0.   0.   0.   0.02 0.12 0.5  0.31 0.04 0.   0.  ]     [0.07 0.   0.   0.   0.   0.   0.26 0.33 0.2  0.15]     [0.73 0.12 0.   0.   0.   0.   0.   0.02 0.02 0.12]]```  B. Algorithme de décodage, ViterbiRetrouver **la sequence d'état caché la plus probable** ayant permis de générer la séquence d'observation.Rappels sur l'algorithme Viterbi (1967):- Il sert à estimer la séquence d'états la plus probable étant donnés les observations et le modèle.- Il peut servir à approximer la probabilité de la séquence d'observation étant donné le modèle. 1\. Initialisation (avec les indices à 0 en python): $$\begin{array}{ccccccccc} \delta_{0} (i) &=& \log \pi_{i} +\log b_{i} (x_{1}) \\ \Psi_{0}(i) &=& -1 \mbox{ Note: -1 car non utilisé normalement} \end{array}$$2\. Récursion: $$ \begin{array}{ccccccccc} \delta_{t} (j) &=&\displaystyle \left[\max_{i} \delta_{t-1}(i) + \log a_{ij}\right] + \log b_{j}(x_{t}) \\ \Psi_{t}(j) &=&\displaystyle \arg\max_{i\in [1,\ N]} \delta_{t-1} (i) + \log a_{ij} \end{array}$$3\. Terminaison (indices à {$T-1$} en python) $$ S^{\star} = \max_{i} \delta_{T-1}(i)$$4\. Chemin $$\begin{array}{ccccccccc} s_{T-1}^{\star} & = &\displaystyle \arg\max_{i} \delta_{T-1}(i) \\ s_{t}^{\star} & = & \displaystyle \Psi_{t+1}(s_{t+1}^{\star}) \end{array}$$L'estimation de $\log p(x_0^{T-1} | \lambda)$ est obtenue en cherchant la plus grande probabilité dans la dernière colonne de $\delta$. Donner le code de la méthode `viterbi(x,Pi,A,B):`  B.1 Définition de la fonction viterbi<jupyter_code>def viterbi(x,Pi,A,B,epsilon=1e-5):
    T = len(x)
    N = len(Pi)
    logPI = np.log(Pi + epsilon)
    logA = np.log(A + epsilon)
    logB = np.log(B + epsilon)
    logDelta = np.zeros((N,T))
    psi = np.zeros((N,T))
    S = np.zeros(T) # les états à retourner
    
    # 1) Initialisation
    for i in range(K):
        logDelta = logPI + logB[:,i]
    psi[:,0]= -1
        
    # 2) Recursion
    for n in range(1, K):
        for i in range(N):
            tmp_sum = logA[:, i] + psi[:, n-1]  # somme temporaire pour reutiliser dans lignes suivantes
            logDelta[i, n-1] = np.argmax(tmp_sum)
            psi[i, n] = np.max(tmp_sum) + logB[i, int(x[n])]

    # 3&4) Terminaison et Chemin
    S = np.zeros(K)
    S[-1] = np.argmax(psi[:, -1])
    for i in range(K-2, 0, -1):
        S[i] = logDelta[int(S[i+1]), i]
        
    # 5) Calcul de la log vraissemblance
    print(logDelta)
    logL = np.max(logDelta[:,-1])
        
    return S, 0 # états + approximation de la vraimsemblance des obs par rapport au modèle
def viterbi(allx,Pi,A,B):
    etat = A.shape[0]    
    taille= len(allx)  
    epsilon = 10e-5   # valeur petite pour empecher log(0) error
    A_log = np.log(A + epsilon)
    Pi_log = np.log(Pi + epsilon)
    B_log = np.log(B + epsilon)

    # 1) Initialisation
    psi = np.zeros((len(A), len(allx))) # A = N
    psi[:,0]= -1
    delta = np.zeros((len(A), len(allx)))  # initialisation en dimension mais pas en contenu !
    
    # 2) Recursion
    for n in range(1, taille):
        for i in range(etat):
            temp_sum = A_log[:, i] + psi[:, n-1]  # somme temporaire pour reutiliser dans lignes suivantes
            psi[i, n] = np.max(temp_sum) + B_log[i, allx[n]]
            delta[i, n-1] = np.argmax(temp_sum)

    # 3&4) Terminaison et Chemin
    S_opt = np.zeros(taille)
    S_opt[-1] = np.argmax(psi[:, -1])
    for n in range(taille-2, 0, -1):
        S_opt[n] = delta[int(S_opt[n+1]), n]

    return S_opt
q, p = viterbi(Xd[0], Pi, A, B) # attention à bien donner la version discrétisée des observations

print("Valeurs:", Xd[0])
print("Decodage:", q)
print("LogL:", p)<jupyter_output><empty_output><jupyter_text>Check:``` décodage[1. 9. 8. 8. 8. 8. 8. 9. 3. 4. 5. 6. 6. 6. 7. 7. 8. 9. 0. 0. 0. 1. 1.]  [0. 0. 0. 0. 0. 0. 0. 0. 1. 2. 2. 2. 2. 2. 3. 3. 3. 3. 4. 4. 4. 4. 4.]  log vraisemblance-38.09356554559258``` B2. Utilisation de l'algorithme de Viterbi pour la classificationNous aurions normalement du utiliser l'algorithme forward pour estimer la vraisemblance... Mais l'algorithme de Viterbi est plus versatile et il nous permettra de mettre en oeuvre des stratégies d'apprentissage dans les questions suivantes.Nous allons donc utiliser la vraisemblance du chemin le plus probable (dans l'espace des états) comme une approximation de la vraisemblance de l'ensemble des chemins.1. Calculer un modèle par lettre1. Estimer la vraisemblance des lettres dans chaque modèle1. Donner un score de classification<jupyter_code># calcul des 26 modèles
models = 

# affectation des signaux = calcul de vraisemblance[n_sig x n_modèles]
all_proba = 
ypred = np.array(all_proba).argmax(1) # à modifier si votre all_proba est transposé

traduction = {ch:i for i,ch in enumerate(np.unique(Y))}
Ynum = np.array([traduction[y] for y in Y])

tx_bonne_classif = np.where(ypred == Ynum,1.,0).mean()
print("Taux de bonne classification" , tx_bonne_classif)<jupyter_output><empty_output><jupyter_text>Check : 85.07% de bonne classification  B3. Séparation apprentissage / testL'expérience précédente n'est pas très satisfaisante, l'évalution des performances n'étant pas fiable.Nous allons donc séparer les données en un jeu d'apprentissage et un jeu de test.<jupyter_code># TOUT LE CODE EST FOURNI ICI

# séparation stratifiée par classe (non randomisée)
def separation_app_test(Y, pc_app):
    ind_app  = []
    ind_test = []
    for cl in np.unique(Y):
        index = np.where(Y == cl)[0]
        ind_app += index[:int(len(index)*pc_app)].tolist()
        ind_test += index[int(len(index)*pc_app):].tolist()
    return ind_app, ind_test

pc_app = 0.8 # pourcentage de points en apprentissage
ind_app, ind_test = separation_app_test(Y, pc_app)
# separation des données
Xd_app, Xd_test = Xd[ind_app], Xd[ind_test]
Y_app, Y_test   = Y[ind_app], Y[ind_test]
Yn_app, Yn_test = Ynum[ind_app], Ynum[ind_test]
Q_app, Q_test   = Q[ind_app], Q[ind_test]
# apprentissage des modèles sur les données d'apprentissage seulement

models = 

# calcul des affectations
all_proba_app = 
all_proba_test = 

# calcul des performances
ypred_app  = np.array(all_proba_app).argmax(1) # a verifier en fonction de vos amtrice all_proba
ypred_test = np.array(all_proba_test).argmax(1)
tx_bonne_classif_app = np.where(ypred_app == Yn_app,1.,0).mean()
tx_bonne_classif_test = np.where(ypred_test == Yn_test,1.,0).mean()
print("Taux de bonne classification" , tx_bonne_classif_app, tx_bonne_classif_test)<jupyter_output><empty_output><jupyter_text>Check :Taux de bonne classification 0.8509615384615384 0.7666666666666667  B4. Matrice de confusionAfin de mieux comprendre les erreurs, tracer la matrice de confusion associée aux lettre<jupyter_code>def mat_conf(Y, Yhat):
    nCl = len(np.unique(Y)) # on espere que le vecteur des Y contient toutes les valeurs entre 0 et nCl-1
    conf = np.zeros((nCl,nCl))
    for y,yh in zip(Y,Yhat):
        conf[y, yh] += 1
    return conf

# appel à la fonction<jupyter_output><empty_output><jupyter_text>B5. OPT. Trouver les paramètres de discrétisation et de nombre d'états optimaux vis à vis des données de testCette procédrure s'appelle un grid-search: proposer des valeurs puis tester de manière exhaustive. Cette procédure est plutot simple à implémenter... mais longue en temps de calcul.**NE PAS FAIRE EN TP****Note:** à force de faire des évaluations sur les données de test, on introduit un risque de sur-apprentissage de ces données.  C. Baum-Welch simplifiéEn utilisant la procédure de Baum-Welch simplifiée vue en TD et rappelée ci-dessous, proposer un code pour apprendre (ou plutot raffiner) les modèles correspondant aux 26 lettres de l'alphabet.**Baum-Welch simplifié:**1. Initialisation des états cachés arbitraire (eg méthode gauche-droite)1. Tant que ''critère de convergence'' non atteint    1. Apprentissage des modèles $\lambda_{lettre}=\{\Pi, A, B\}$    1. Estimation des états cachés par ViterbiLe critère de convergence sera la vraisemblance.* A chaque itération $k$ et pour toutes les lettres $lettre$, calculer pour l'ensemble des séquences d'observation : $$\log\mathcal L^k = \sum_{lettre}\sum_i \log p(\mathbf x_i^{lettre} | \lambda_{lettre}^k)$$* Lorsque la vraisemblance n'évolue plus (ie $\frac{\log\mathcal L^k - \log\mathcal L^{k+1}}{\log\mathcal L^k} < 1e-4$), sortir de la boucle de mise à jour.* Donner l'implémentation de la méthode d'apprentissage* Tracer la courbe de l'évolution de la vraisemblance au cours des itérations C1. Donner le code de la fonction baumwelch<jupyter_code># nIter : nombre d'iteration max (pour éviter les boucles infinies)
# Xd : les données d'observations discrétisées
# Yn : les étiquettes au format numérique
# models: une initialisation maline des modèles à apprendre
#     - RAPPEL: l'initialisation est critique !

def baumwelch(nIter, Xd, Yn, models):
    N = len(models[0][0])
    K = models[0][2].shape[1] # la seconde dimension de B
    print(N,K)
    # Baum Welch (simplifie)
    nCl = len(np.unique(Y))
    nIter = 10
    # A compléter
    
    return models

models = baumwelch(30, Xd_app, Yn_app, models)<jupyter_output><empty_output><jupyter_text>C2. Evaluation des performancesLes performances sont-elles meilleures après optimisation du modèle?Faudrait-il refaire une procédure de grid-search basée sur l'algorithme Baum-Welch?<jupyter_code># calcul des affectations en apprentissage et en test
all_proba_app = 
all_proba_test = 

# calcul des perfromances
ypred_app  = np.array(all_proba_app).argmax(1)
ypred_test = np.array(all_proba_test).argmax(1)
tx_bonne_classif_app = np.where(ypred_app == Yn_app,1.,0).mean()
tx_bonne_classif_test = np.where(ypred_test == Yn_test,1.,0).mean()
print("Taux de bonne classification" , tx_bonne_classif_app, tx_bonne_classif_test)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13807837870833842, 'fraction_numerical': 0.06635331461360029, 'mean_word_length': 3.2564302416212003, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '383044', 'license': 'MIT', 'n_tokens_mistral': 7177, 'n_tokens_neox': 6552, 'n_words': 1927, 'path': 'S1/MAPSI/TME/TME7/enonce2/TME7-HMM-lettres-DurandKordon.ipynb', 'repo_name': 'hanzopgp/MasterBioInformaticsPTs', 'size': '16382'}