////////////////////////////////////////////////////////////////////////////////////
////        Equations de Dyson                                                  ////
////       (champ moyen répulsif, Algorithme MALA (schéma d'euler + étape MCMC) ////
////////////////////////////////////////////////////////////////////////////////////


clf();
stacksize("max");
T=1;
n=100; //Pas de temps du schéma d'Euler h=T/n
time_mesh=(T/n)*[0:n]';
N=100;  //Nombre de particules 10<N<100 (pour N=100 régler la condition initiale avec epsilon)
 
nmax=10000; //(nmax 10 000 car boucle lente)  Exemple : n=100, 10<N<100, nmax=10000


///////////////////////////////////////////////////////////////////////////////////////////
//Condition initiale = N points entre - eps (N sqrt(h)/2) et  - eps (N sqrt(h)/4)   ///////
///////////////////////////////////////////////////////////////////////////////////////////

eps=.5;  // eps=1,1.5, 2
X=eps*[1:N/2]/(sqrt(n/T))
X=[gsort(-X,'lc','i'),X];

pX=[];
VX=[];
for i=1:N
U=X-X(i);
V=((U>0)|(U<0));
pX=[pX,abs(prod(U(V)))]; //Déterminant de Vandermonde
VX=[VX,(1/N)*sum(U(V)^(-1))]; //Dérive initiale du processus de Dyson
end



Z=X;

//s=size(Z); 


temps_reel=1;

for t=1:nmax

//////////////////////////////////////////////////////////////////////////////////////////////////
////// COMPLETER les lignes rouges :  EVOLUTION DU MODELE DE DYSON APRES L'ETAPE MCMC         ////
//////////////////////////////////////////////////////////////////////////////////////////////////
        
Y=X+	  ;


//Début de la transition de MCMC (Algorithme MALA)

Ysort=gsort(Y,'lc','i'); 
test=sum(double(((Y-Ysort)==0))); // Vérfification d'appartenance à la chambre de Weyl

if test==N then

// Calcul des déterminants de Vandermonde et des dérives de Dyson pour l'état proposé
pY=[];
VY=[];
for i=1:N
U=Y-Y(i);
V=((U>0)|(U<0));
pY=[pY,abs(prod(U(V)))];
VY=[VY,(1/N)*sum(U(V)^(-1))];
end


// Paramètre d'acceptation rejet de l'algorithme de Metropolis-Hasting )
r=exp(log(prod(pY./pX))+(N/(4*(temps_reel+(T/n))))*sum(X^2-Y^2) +(N/2)*sum(((X-Y).*(VX+VY)+T/(2*n)*(VX^2-VY^2))));
a=min(1,r);

// Etape d'acceptation rejet
if  grand(1,1,'def')<=a then
X=Y;
VX=VY;
pX=pY;
Z=[Z;Y];
temps_reel=temps_reel+(T/n)
end




end
end


// Affichage des trajectoires


clf(0);xset("window",0);show_window();
  
s=size(Z); 
    for i=1:N
        plot2d((T/n)*[1:s(1)]',Z(:,i),i);
    end