//////////////////////////////////////////////////////////////////////////////
//   Simulation de Marches aléatoires simples sur Z^2 sans intersection
//   par des arbres généalogiques de systèmes de particules en interaction.
//   Estimation des constantes de connectivité.
//
//
/////////////////////////////////////////////////////////////////////////
clear
stacksize(90000000)
T=100;
U=[1 0 -1 0 ; 0 1 0 -1];         // Vecteurs unité sur Z^2
I=floor(1+4*grand(1,T,'def'));
X=zeros(2,1);
//Simulation d'une marche aléatoire simple
for t=1:T
X=[X,X(:,t)+U(:,I(t))];
end
xset("window",1);
xbasc();
plot2d(X(1,:),X(2,:),5,rect=[-sqrt(2*T),-sqrt(2*T),sqrt(2*T),sqrt(2*T)]);
 xtitle(['Marche aleatoire simple sur Z^2, Horizon temporel :' +string(T)],'','')




N=100;                           //Nombre de particules 
J=floor(1+4*grand(T,N,'def'));   //Choix aléatoires des vecteurs directeurs
UT=[];
for t=1:T
 UJ=[];
 for i=1:N
 UJ=[UJ;U(:,J(t,i))];
 end
UT=[UT,UJ];
end

XI=zeros(2*N,1);                  //Initialisation des marches en l'origine
for i=1:N
YI(1:2,1,i)=[XI(2*i-1,1);XI(2*i,1)];
end


//QQ variables initiales
G=zeros(N,T);
Proba=1;
LS=[];
A=1;   //Pour tester l'arret de l'algorithme lorsque toutes ont intersecte
for t=2:T

   //Explorations aléatoires élémentaires des N marches  
   for i=1:N
   YI(:,1:t,i)=[YI(:,1:t-1,i),YI(:,t-1,i)+UT(2*i-1:2*i,t)];
   G(i,1)=[0];
   P(i)=1;
     for s=2:t // Calcul des distances L1 entre les états courants 
               // et les points visités par le passé    
     G(i,1:s)=[G(i,1:s-1),abs(YI(1,t,i)-YI(1,s-1,i))+abs(YI(2,t,i)-YI(2,s-1,i))];
     end
   G(i,1:t)=[G(i,2:t),0];
   //P(i)=prod(G(i,1:t-1));
    Q=-sort(-G(i,1:t-1));
    P(i)=Q(1);
  end

     Inter=find(P==0);   //Les indices des configurations présentant une intersection
     LInter=length(Inter); 
     NInter=find(P>0);   //Les indiuces des configurations sans intersection
     LNInter=length(NInter);
     if Inter==N then A=0
     end            
     Proba=[Proba,Proba(t-1)*(LNInter/N)]; //Estimation de la proba de  non intersection
                                           //de la marche aléatoire simple après t itérations
     LS=[LS,log(4^t*Proba(t))/t];          //Estimation du logarithme de la constante 
                                           //de connectivité
     //Sélection des configuration sans intersection  
     ZI=YI(:,1:t,NInter);
     V=floor(1+LNInter*grand(1,LInter,'def'));
     if Inter<>[] then
     YI(:,1:t,Inter)=ZI(:,1:t,V);
     end
     
D=max(abs(YI(:,1:t,:))); //Réglage des dimension des fenêtres


  xset("window",2);
  xbasc();
  for i=1:N
  plot2d(YI(1,1:t,i),YI(2,1:t,i),i,rect=[-D-1,-D-1,D+1,D+1]);
  end
  xtitle(['Evolution de l arbre genealogique, horizon temporel :' +string(t);
           'Nombre de particules :' +string(N)],'','')



  xset("window",3);
  xbasc();
  plot2d(YI(1,1:t,1),YI(2,1:t,1),1,rect=[-D-1,-D-1,D+1,D+1]);
   xtitle(['Evolution de la premiere ligne ancestrale, horizon temporel :' +string(t);
           'Nombre de particules :' +string(N)],'','')

end

xset("window",4);
xbasc();
for i=1:N
plot2d(YI(1,1:T,i),YI(2,1:T,i),i,rect=[-D-1,-D-1,D+1,D+1]);
end
xtitle(['Arbre genealogique, horizon temporel :' +string(T);
           'Nombre de particules :' +string(N)],'','')


xset("window",5);
xbasc();
plot2d(YI(1,1:T,1),YI(2,1:T,1),1,rect=[-D-1,-D-1,D+1,D+1]);
 xtitle(['Premiere ligne ancestrale, horizon temporel :' +string(T);
           'Nombre de particules :' +string(N)],'','')



xset("window",6);
xbasc();
plot2d(2:T,Proba(2:T),1);
xtitle(['Estimation de la probabilite de non intersection, nombre de particules :' +string(N)],"axe temporel","probabilite")




xset("window",7);
xbasc();
plot2d(1:T-1,log(2)*ones(1:T-1),3);
plot2d(1:T-1,LS(1:T-1),2);
plot2d(1:T-1,log(3)*ones(1:T-1),3);
xtitle(['Estimation du log de la constante de connectivite, nombre de particules :' +string(N)],'axe du temps','')