/////////////////////////////////////////////////////////////////////////////////
// Simulation par des arbres généalogiques d'une marche aléatoire sur Z 
// absorbee au bords d'un segment [-L,L], conditionnellement à la non absorption. 
// Estimation particulaire de la probabilite de survie.  
// Estimation de la valeur propre principale et du vecteur propre associé (cas réversible)
// d'une matrice d'absorption positive.
//
//
/////////////////////////////////////////////////////////////////////////////////
//QQ Valeurs :
//(T,N,p,L)=(100,5,0.5,5) bonne vision de l'arbre
//(T,N,p,L)=(100,10,0.5,7) bonne vision de l'arbre
//(T,N,p,L)=(100,10000,0.5,4 ou 7) bonne vision du vecteur propre+estimee v.p.
//(T,N,p,L)=(1000,1000,0.5,7) bonne vision du vecteur propre+estimee v.p.

//(T,N,p,L)=(1000,1000,0.5,5) estimee de la v.p
//(T,N,p,L)=(1000,1000,0.5,5) estimee de la v.p
//(T,N,p,L)=(1000,1000,0.7,7) estimee de la v.p
//(T,N,p,L)=(400,10000,0.7,3) estimee de la v.p
//(T,N,p,L)=(1000,1000,0.7,10) estimee de la v.p
//(T,N,p,L)=(1000,1000,0.7,7) estimee de la v.p
/////////////////////////////////////////////////////////////////////////////////

clear;
T=100;                                  //Horizon temporel
N=10000;                                   //Nombre de particules
p=0.5;                                   //Proba( X(n)=x+1 ! X(n-1)=x)
L=4;                                     //Niveau absorbant
proba=1;
//x=zeros(N,T);
//t=[1:T];
u=zeros(N,T);
x=floor(1+(2*L+1)*rand(N,T,'def')-L);





// Description des matrices d'absorption
d=2*L+1;
P=p*ones(1,d-1);
Dp=diag(P);
M1=[zeros(d,1),[Dp;zeros(1,d-1)]];
Dq=diag((1-p)*ones(1,d-1));
M2=[zeros(1,d);[Dq,zeros(d-1,1)]];
M=M1+M2;
g=[0,ones(1,d-2),0];
G=diag(g);
Q=G*M;

//Calcul du log de la valeur propre maximale
vp=sort(spec(Q));
lvp=log(vp(1));


//Estimation du log de la valeur propre principale
//par des itérations de la matrice d'absorption
I=diag(ones(1,d));
lvpn=[];
I=diag(ones(1,d));
for s=1:T
I=I*Q;
I1=sort(I*ones(d,1));
lvpn=[lvpn,log(I1(1))/s];
end





lvpp=[];
hv=zeros(d,T);


for s=2:T

    //Evolution des particules
    u(:,s)=bool2s(rand(N,1,'uniform')<p);
    x(:,s)=x(:,s-1)+2*u(:,s)-1;

//Calcul dees taux de réussites=proportions de non absorption   
w=find(abs(x(:,s))>=L);
l=length(w);
v=find(abs(x(:,s))<L);
q=length(v);
proba=proba*(q/N);

lvpp=[lvpp,log(proba)/s];  //Estimation particulaire du log de la valeur propre principale


//Mise à jour/sélection des particules non absorbées
     y=x(v,1:s);
     U=floor(1+q*grand(1,l,'def')); 
     if w<>[] then
     x(w,1:s)=y(U,1:s);
     end

//Mesure d'occupation spatio-temporelle (pour l'estimation du vecteur propre
//dans le cas réversible)
if p==0.5 then
for m=-L:L
c=find(x(:,s)==m);
hv(m+L+1,1:s)=[hv(m+L+1,1:s-1),length(c)/N];
end
end



//Evolution de l'arbre généalogique
    // xset("window",2);
     //xbasc();
     //plot2d(1:s,x(:,1:s)',rect=[0,-L,T,L])
     //xtitle("Evolution de l arbre genealogique","axe temporel","Z")

     
          
end

//Arbre généalogique à l'instant terminal
//xset("window",2);
//xbasc();
//plot2d(1:T,x(:,1:T)',rect=[0,-L,T,L])
//xtitle("Arbre genealogique terminal","axe temporel","Z")


//Repartition terminale des particules (=vecteur propre principal 
//dans le cas réversible, p=1/2) 
if p==0.5 then
 xset("window",3);
 xbasc();

hT=[];
for m=-L:L
hT=[hT,sum(hv(m+L+1,:))/T];
end

plot2d(-L:L,hv(:,T)',style=2);
plot2d(-L:L,hT',style=3);
if p==1/2 then [evals,X]=spec(Q);
VP=evals(:,1).*sign(evals(:,1));
VPnorm=VP/sum(VP);
plot2d(-L:L,VPnorm',style=5);
end
xtitle("Estimation du vecteur propre principal (theorique=rouge, particulaire (instant terminal=bleu, spatio-temporelle=vert))","Z","")
end



xset("window",1);
xbasc();
xselect();
lvpo=lvp*ones(1,T);
plot2d([1:T],lvpn,3)
plot2d([1:T],lvpo,5)
plot2d([2:T],lvpp,2)
xtitle("Logarithmes des valeurs propres (theorique=rouge, matrices iterees=vert, particulaire=bleu)","axe temporel","")


xset("window",4);
xbasc();
xselect();
plot2d([2:T],lvpn(2:T)-lvpp,1)
xtitle("Evolution des erreurs d estimation des v.p.","axe temporel","")


lvp
lpvT=lvpn(T)
lvpa=log(proba)/T