import numpy as np
import matplotlib.pyplot as plt

M=10000 # taille echantillon Monte Carlo
I=10 # nombre de types de sites
s=20 # population totale
it=10 # nombre initial d'iterations


lam=10*np.random.random(I) # 1+10*rand(1,I,'u') parametres des lois de Poisson

#initialisation des nombres de sites de chaque type occupes



X=np.zeros((M,I)) #vecteur des nombres de sites de chaque type occupes

for m in range(M):
    Y=np.floor(I*np.random.random(s)); #Tirage uniforme du type de site occupe par chacune des s particules
    for l in range(s):
        X[m,int(Y[l])]=X[m,int(Y[l])]+1 #comptage des nombres de sites de chaque type occupesY=1+floor(I*rand(M,s,'u')); //Tirage uniforme du type de site occupe par
                            
#X[:,0]=20*np.ones(M) initialisation deterministe

m=np.amax(lam) #maximum du vecteur lambda
l=np.argmax(lam)  #indice de probabilite maximale dans la loi multinomiale
p=lam[l]/np.sum(lam) #probabilite correspondante
q=1-p
pr=q**s
pro=np.ones(s+1)
for i in range(s+1):
      pro[i]=pr
      pr=pr*p*(s-i)/(q*(i+1))
#boucle sur le nombre d'iterations
n=it
for nbit in range(3):
  for k in range(it):
    #tirage de M couples (i,j)
    i=np.sum(np.cumsum(X,axis=1)<s*np.random.random((M,1))*np.ones((1,I)),axis=1)
    j=np.floor((I-1)*np.random.random(M)).astype(int)
    j=j+(j>=i)
    #mise a jour des M etats 
    inc=(np.random.random(M)<lam[j]/lam[i])
    for m in range(M):
        X[m,i[m]]=X[m,i[m]]-inc[m]
        X[m,j[m]]=X[m,j[m]]+inc[m]
  

  compt=np.zeros(s+1);
  #comptage dans compt(x) du nombre de simulations m pour lesquelles
  #X(m,l)=x
  
  

  plt.figure()
  plt.bar(range(s+1), compt/M, width = 0.1, color='r',label="frequences empiriques n="+ str(n) )
  #Trace en rouge le diagramme en baton des moyennes empiriques des
  #nombres de partcules qui occupent les sites de probabilite maximale
  
  plt.bar(np.arange(s+1)+0.2,pro, width = 0.1, color='b',label="loi binomiale limite")
  plt.legend(loc="best")
  it=n*3 #multiplication par 3 du nombre d'iterations supplementaires
  n=n*4 #multiplication par 4 de l'horizon