/////////////////////////////////////////////////////////////////
// 
// Proba d'atteinte d'un niveau supérieur avant d'etre absorbee en un niveau inferieur
// (Marche aléatoire X sur Z, Proba(X(n)=x+1! X(n-1)=x)=p)
// Simulation/Estimation par des arbres généalogiques.
//
//////////////////////////////////////////////////////////////////////////////////////////
//QQ Valeurs : T=100 (ou 200 ou 300), N=50 ou 100 ou 150, Niv=de 7 a 10, p= de 0.1 a 0.9 
//par pas de Dalpha=alpha0=1, 2,  et etat absorbant L=-2                                      
//(T,N,Niv)=(20,5,5),(200,150,10 par pas de 2)
/////////////////////////////////////////////////////////////////////////////////////////

clear
stacksize(90000000)
T=100; //longueur maximale des excursions
N=150;  //Nombre de particules
Niv=4; //Nombre de niveaux
p=0.1;   //Proba(X(n)=x+1! X(n-1)=x)=p
L=-2;    //Niveau obstacle

x0=0;     //Initialisation des particules
alpha0=1;  //Niveau initial

x(:,:,1)=x0*ones(N,T+1); //Initialisation des 1res excursions
Texc(:,:,1)=ones(N,T+1);  //Initialisation des horloges des 1res excursions.

alpha=alpha0*ones(1,T+1); 


Tmax=ones(1,Niv);


//Description des niveaux

Dalpha=2;
for n=1:Niv
alpha(n+1)=alpha(n)+Dalpha;
end

//Initialisations de proba. de passages/non absorption

proba=ones(1,Niv+1);


for n=1:Niv // Iterations des excursions entre niveaux d'indices n

lv=0; //Nombre de particules aborbees
lw=0; //Nombre de particules ayant franchi le niveau
s=1;  //Initialisation du compteur

   while lw+lv<N
    
             //Procedure d'arret des particules absorbees    
             w=find(x(:,s,n)<= L); 
             lw=length(w);
             if w<>[] then x(w,s+1:T,n)=x(w,s,n)*ones(1,T-s); 
                           Texc(w,s+1:T,n)=Texc(w,s,n)*ones(1,T-s);
             end
             
             //Procedure d'arret des particules ayant franchi le niveau alpha(n)
             v=find(x(:,s,n)>= alpha(n)); 
             lv=length(v);
             if v<>[] then x(v,s+1:T,n)=x(v,s,n)*ones(1,T-s); 
                           Texc(v,s+1:T,n)=Texc(v,s,n)*ones(1,T-s);
             end             
             
      //Evolution de particules restantes
      if lw+lv<N then    
                     for i=1:N
                       if x(i,s,n)<alpha(n) then if x(i,s,n)>L then x(i,s+1,n)=x(i,s,n)+2*bool2s(rand(1,1,'uniform')<p)-1;   
                                                                Texc(i,s+1,n)=Texc(i,s,n)+1;
                                                           end                
                                           end
                     end           
      end

 
 
     xset("window",1);
     xbasc();
     plot2d(1:s,x(:,1:s,n)')
     plot2d(1:s,alpha(n)*ones(1,s))
     plot2d(1:s,L*ones(1,s))
     xtitle("Exploration du niveau : "+string(alpha(n)),"","")


     
    s=s+1; //Increment du compteur
    end 
 
Tmaxm=sort(Texc(:,:,n));
Tmax(n)=Tmaxm(1,1);    //Long. excursion la plus longue
 
 
//Boucle de calcul de proba. de succes de passage des niveaux/non absorption 
proba(n+1)=proba(n)*(length(v)/N);
      
      
      //Mise en memoire et selection des excursions ayant passe le niveau n.
      y=x(v,1:T+1,n); TTexc=Texc(v,1:T+1,n);
      U=floor(1+length(v)*grand(1,length(w),'def'));
      if w<>[] then x(w,1:T+1,n)=y(U,1:T+1);
                    Texc(w,1:T+1,n)=TTexc(U,1:T+1);
      end    

       

     xset("window",2);
     xbasc();
     plot2d(1:Tmax(n),x(:,1:Tmax(n),n)')    
     plot2d(1:Tmax(n),alpha(n)*ones(Tmax(n),1),5)
     xtitle("Excursions selectionnees au niveau : "+string(alpha(n)),"","")


//Re-initialisation des excursions et des horloges   
   
x(:,:,n+1)=alpha(n)*ones(N,T+1);
Texc(:,:,n+1)=ones(N,T+1);

   
end    


//Calcul des proba. theoriques de passage de x0 au niveau n sans toucher L   
probth=ones(1,Niv+1);
q=1-p;

if p<>.5 then
             for n=1:Niv
             probth(n+1)=((q/p)^(x0)-(q/p)^L)/((q/p)^(alpha(n))-(q/p)^(L));
             end
end

if p==0.5  then         
             for n=1:Niv probth(n+1)=(x0-L)/(alpha(n)-L);
             end      
end


     xset("window",4);
     xbasc();
     plot2d(1:Niv+1,probth(1:Niv+1),5)    
     plot2d(1:Niv+1,proba(1:Niv+1),2)   
     xtitle("Proba. d atteinte avant absorption (theorique en rouge), nombre de part. : "+string(N),"Nb. de niveaux","")


     
     
     xset("window",5);
     xbasc();
     plot2d(1:Niv,Tmax(1:Niv),5)    
     xtitle("Durees maximales des excursions entre niveaux","niveaux","")

       
       
     //Calcul des durees de chaque excursions entre niveau
     TEXC=zeros(N,Niv);
     for i=1:N 
     for n=1:Niv
     F=find(Texc(i,:,n)==[1:T+1]);
     Fm=sort(F);
     TEXC(i,n)=Fm(1);
     end
     end
     
     
       xset("window",6);
     xbasc();
    for i=1:N
     Exc=x(i,1:TEXC(i,1),1);
     for n=2:Niv
     Exc=[Exc,x(i,2:TEXC(i,n),n)];
     end
     plot2d(Exc);
     xtitle("Arbre genealogique complet des particules selectionnees, nb de part. : "+string(N),"","")
     end
     
     
  
     xset("window",7);
    xbasc();HExc=[];
    for n=1:Niv
      for i=1:N
      HExc=[HExc;x(i,1:TEXC(i,n),n)'];
      end
     //MP=[MP;Exc]; 
     histplot(10,HExc',n);
     xtitle("Distribution des mesures d occupation entre niveaux (n=1 noir, n=2 bleu, n=3 vert,...), Nb part.: "+string(N),"niveaux","")
    end
          
     
     xset("window",8);
     xbasc();
     TMOY=[];
     for n=1:Niv
     TMOY=[TMOY,sum(TEXC(:,n))/N];  
     end  
     plot2d(1:Niv,TMOY,5)    
     xtitle("Durees moyennes des excursions entre niveaux","niveaux","")