"""
Simulation d'une chaine de Markov
"""
import numpy as np
import matplotlib.pyplot as plt

  

def Markov(A,x):
    """
    Markov(A,x): effectue un tirage aléatoire et fait évoluer la chaine 
    dont l'état initial est x, retourne le nouvel état 
    """
    p=A[:,x] #extrait la colonne x de A
    choix=np.arange(A.shape[0]) #une liste des entiers où faire le choix
    return np.random.choice(choix,p=p) #retourne y de choix avec proba p[y]

#Version sans random.choice, uniquement avec random.rand()
#(cf TP simulation aléatoire)
def Markov(A,x):
    """
    Markov(A,x): effectue un tirage aléatoire et fait évoluer la chaine 
    dont l'état initial est x, retourne le nouvel état 
    """
    p=A[:,x] #extrait la colonne x de A
    u=np.random.rand()
    q=0
    for k in range(A.shape[0]):
        q+=p[k]
        if u <= q:
            return k

    
#fonction principale
def simulation_Markov(N,A,x0,NS=1000) :
    """    
    simulation_Markov(N,A,x0) :  fait NS simulations de la chaine de markov 
        commençant à l'état x0  avec matrice de transition A
        affiche l'histogramme empirique des états pris
        affiche aussi l'histogramme théorique obtenu par multiplication matricielle
        donne la différence maximale des probabilités de présence entre la simulation et le calcul
        retourne les probabilités de présence calculées
    """
    #calculs théoriques
    P=np.zeros((A.shape[0],)) # fabrique un vecteur de zeros autant de lignes que A
    P[x0]=1.0 # place l'état initial
    n=N # on calcule la puissance N-ieme
    while n>0 :
        n=n-1
        P=np.dot(A,P)
    P_calcule=P




    #liste pour compter le nombre de passages dans chaque état, contient autant d'éléments que le tableau A a de lignes
    P=np.zeros((A.shape[0],))
    s=NS #nbre de simulations
    while s>0:
        s=s-1
        x=x0 #on initialise la chaine au départ donné en paramètre
        n=N # on fait maintenant tourner la chaine pendant N tours
        while n>0 :
            n=n-1
            x=Markov(A,x)

        P[x]+=1 # on ajoute 1 à la où on est arrivés
    
    P_simule=P/NS # fréquences de passage
    
    plt.figure()
    plt.xlim(-0.25,A.shape[0]-0.75)
    plt.ylim(-0.05,1.05)
    plt.xticks(range(A.shape[0]))
    plt.plot(P_simule,'rx') # simulé en rouge avec des croix
    plt.plot(P_calcule,'bo') # calculé en bleu avec des ronds

    
    print("MARKOV : Chaine donnée par A=",A," x0=",x0,". N=",N,"étapes, NS=",NS," simulations")
    print("la différence maximale de densité entre la simulée et la théorique est de ", np.max(np.abs(P_calcule-P_simule)))
    return P_simule    
    
#programme principal

#initialise la chaine, met en place la matrice de transition A_ij = P(X_1=i|X_0=j)
#partant de 0
A=np.array([
[0.5, 0.2, 0.3],
[0.3, 0.4, 0.3],
[0.2, 0.4, 0.4],
])
x0=0 # état initial, les états sont numéroté de 0 à A.shape[0]-1

d=simulation_Markov(10,A,x0)
plt.savefig("simulation-markov-01.pdf",format='pdf')


A=np.array([
[0.5, 0.2, 0.3],
[0.1, 0.4, 0.3],
[0.4, 0.4, 0.4],
])
x0=0 # état initial, les états sont numéroté de 0 à A.shape[0]-1

d=simulation_Markov(10,A,x0)
plt.savefig("simulation-markov-02.pdf",format='pdf')