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
kappa=np.zeros(M) #somme des valeurs de X(l) au cours des iterations

#boucle sur le nombre d'iterations
n=it
for nbit in range(4):
  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]
    #######################
    # mettre à jour kappa #
    #######################
    
  plt.subplot(2,2,nbit+1)

  ########################################################################## 
  #calculer moy=moyenne empirique de kappa/n et Z=moy+sqrt(n)(kappa/n-moy) #
  ##########################################################################
  


  plt.hist(Z,density = True,bins=20,label="n="+ str(n))
  v=np.var(Z)
  Z=np.sort(Z)
  plt.plot(Z,np.exp(-(Z-moy)**2/(2*v))/np.sqrt(2*np.pi*v))
  plt.legend(loc="best")
  it=n*3#multiplication par 3 du nombre d'iterations supplementaires
  n=n*4#multiplication par 4 de l'horizon