clf 
// fonctions de base
function S1=tau0(S,i,j)
S1=S
S1(i,j)=0
endfunction

function S1=taum1(S,i,j)
N=size(S)
N=N(1)
S2=zeros(N+2,N+2)
S2(2:(N+1),2:(N+1))=S
if or([S2(i+2,j+1)==1,S2(i,j+1)==1,S2(i+1,j+2)==1,S2(i+1,j)==1]) then
    S1=S
else
    S1=S
    S1(i,j)=-1    
end
endfunction


function S1=taup1(S,i,j)
N=size(S)
N=N(1)
S2=zeros(N+2,N+2)
S2(2:(N+1),2:(N+1))=S
if or([S2(i+2,j+1)==-1,S2(i,j+1)==-1,S2(i+1,j+2)==-1,S2(i+1,j)==-1]) then
    S1=S
else
    S1=S
    S1(i,j)=1    
end
endfunction

\\ Fonction de transition pour la chaine de Markov

function S1=trans(S,q,i,j,u)
if u<(-1/(2*q+1))
    S1=taum1(S,i,j)
elseif u>(1/(2*q+1))
    S1=taup1(S,i,j) 
else
    S1=tau0(S,i,j)
end
endfunction


\\ Une trajectoire de la chaine de Markov

N=10;
S=zeros(N,N);
q=1;

for i=1:100
    i=grand(1,1,'uin',1,N);
    j=grand(1,1,'uin',1,N);
    u=2*rand()-1;
    S=trans(S,q,i,j,u);
end
Matplot(2*(S+1))

\\ Coupling from the past
N=10;
U=ones(N,N);
V=-U
q=1;
n=0;
i=grand(1,1,'uin',1,N);
j=grand(1,1,'uin',1,N);
u=2*rand()-1;

while ~and(U==V)
n=n+1
U=ones(N,N);
V=-U;
i=[i,grand(1,2^n,'uin',1,N)];
j=[j,grand(1,2^n,'uin',1,N)];
u=[u,2*rand(1,2^n)-1];
for k=1:2^n    
U=trans(U,q,i(1,k),j(1,k),u(1,k));
V=trans(V,q,i(1,k),j(1,k),u(1,k));
end    
end

S=U;
Matplot(2*(S+1))

\\ Méthode Monte Carlo Markov Chain pour calculer le nombre moyen de composant en fonction de q
N=10;
moyenne=0
for q=0.1:0.1:2.5
S=zeros(N,N);
for i=1:1000
    i=grand(1,1,'uin',1,N);
    j=grand(1,1,'uin',1,N);
    u=2*rand()-1;
    S=trans(S,q,i,j,u);
end
somme=0;
for i=1:10000
    i=grand(1,1,'uin',1,N);
    j=grand(1,1,'uin',1,N);
    u=2*rand()-1;
    S=trans(S,q,i,j,u);
    somme=somme+sum(abs(S));
end
moyenne=[moyenne,somme/10000];
end
clf
plot(0.1:0.1:2.5,moyenne(1,2:26)/N^2)


