//////////////////////////////////////////////
// fonction generant une marche simple dans Z^2 issue de l'origine
function res=marche2(N)
	e=[0 1 ; 0 -1 ; 1 0 ; -1 0];
	alea=floor(4.*rand(1,N-1))+1;
	aux=[0 0];
		for j=1:N-1
		aux=[aux ; aux($,:)+e(alea(j),:)];
	end
	res=aux;
	//xmin=min(aux(:,1)); xmax=max(aux(:,1)); ymin=min(aux(:,2)) ; ymax=max(aux(:,2));
	//isoview(xmin,xmax,ymin,ymax);
	//plot(aux(:,1),aux(:,2))
	//x=0:0.1:8;
	//plot(sqrt(N)*cos(x),sqrt(N)*sin(x))
endfunction

//////////////////////////////////////////////
// la commande next(i,j,n) specifie les voisins de (i,j) dans la grille n x n
function res=next(i,j,n)
    if ([i j]==[0 0]) then 
        res=[0 1 ; 1 0];
    elseif ([i j]==[0 n-1]) then 
        res=[0 n-2 ; 1 n-1]; 
    elseif ([i j]==[n-1 0]) then 
        res=[n-2 0; n-1 1 ]; 
    elseif ([i j]==[n-1 n-1]) then 
        res=[n-2 n-1; n-1 n-2 ]; 
    elseif ((i=0)&(j<>0))|((i=0)&(j<>n-1))
        res=[i j+1 ; i j-1 ; i+1 j];
    elseif ((i=n-1)&(j<>0))|((i=n-1)&(j<>n-1))
        res=[i j+1 ; i j-1 ; i-1 j];
    elseif ((j=0)&(i<>0))|((j=0)&(i<>n-1))
        res=[i-1 j ; i+1 j ; i j+1];
    elseif ((j=n-1)&(i<>0))|((j=n-1)&(i<>n-1))
        res=[i-1 j ; i+1 j ; i j-1];
    else 
        res=[i-1 j ; i+1 j ; i j-1 ; i j+1];
    end
endfunction


//////////////////////////////////////////////
// la commande whatsnext(i,j,n) choisi uniformement un voisin de (i,j) dans la grille n x n
function res=whatsnext(i,j,n)
    u=next(i,j,n); //on genere la liste de voisins
    k=length(u(:,1));// on compte combien il y en a 
    alea=grand(1,1,'uin',1,k);//on en choisit un au hasard
    res=u(alea,:);// on attribut le resultat
endfunction

//////////////////////////////////////////////
// renvoie une trajectoire de longueur m de la marche simple dans la grille n x n
function res=markovgrille(n,m,a,b)
    z=zeros(m,2);
    z(1,:)=[a b];
    for j=1:1:(m-1)
        z(j+1,:)=whatsnext(z(j,1),z(j,2),n)
    end
    res=z;
    plot2d(z(:,1),z(:,2))
endfunction

////////////////////////////////////////////////
// fonction auxiliere pour la procedure d'effacement des boucles
// Etant donne un couple "valeur" et un matrice "vect" de taille m*2
// cette fonction evalue si la matrice "vect" contient des occurences (des lignes)
// egales à "valeur" en renvoie l'indice de la derniere telle occurence

function res=testegalite(valeur,vect)
	m=length(vect(:,1)); // vect est une matrice m*2
	aux=zeros(m,1);
	for k=1:m
		aux(k,1)=sum(prod(bool2s(vect(k,:)==valeur),2)); // on compte les occurences
	end
	res=sum(bool2s(cumsum(aux)<sum(aux)))+1; // on renvoie l'indice de la derniere occurence
endfunction

// teste si l'element (i,j) est dans V
function res=estdans(i,j,V)
     aux1=(V(:,1)==i).*(V(:,2)==j);
     aux2=find(aux1==1);
     aux3=sum(aux2);
     if (aux3==0) then
         res=0;
     else 
         res=min(aux2);     
     end
endfunction

////////////////////////////////////////////////
// procedure d'effacement de booucle d'une marche plane
// la fonction LE(walk) trace la marche simple walk et sa version a boucles effacees
// et renvoie les coordonnees de la marche a boucles effacees
function res=LE(walk)
	N=length(walk(:,1)); // on calcule la longueur la marche simple
	temp=testegalite([0 0],walk); // on initialise comme dans l'algorithme
	indice=[temp]; // le vecteur indice contiendra les indice de la marche a boucles effacees
	while ((walk(temp,1)<>walk(N,1))|(walk(temp,2)<>walk(N,2)))// on met en place la boucle
	valeur=[walk(temp,1) walk(temp,2)];
	temp=testegalite(valeur,walk);
	indice=[indice ; temp+1];
	temp=temp+1;
    end
	lerw1=walk(indice,1)';
	lerw2=walk(indice,2)';
	// les coordonnees de la marche a boucles effacees sont (lerW1,lerw2)
	l=[length(lerw1)]; // l est la longueur de la nouvelle marche
	xmin=min(walk(:,1)); xmax=max(walk(:,1)); ymin=min(walk(:,2)) ; ymax=max(walk(:,2));
	isoview(xmin,xmax,ymin,ymax);
	plot2d(walk(:,1),walk(:,2)) // on trace la marche initiale
	xset('thickness',2) 
	plot2d(lerw1,lerw2,style=5)// et en gras sa version a boucles effacees
    res=[lerw1 ; lerw2]';
endfunction





// genere et trace une excursion de marche a boucles effacees issue de (a,b) dans la grille n x n jusqu'a ce qu'elle touche T
function res=LERW(a,b,T,n)
    z=[a b];
    tb=0;
    tT=0;
    while (tT==0)
        aux=whatsnext(z($,1),z($,2),n); // prochain point de la marche
        tT=estdans(aux(1),aux(2),T); // on teste si l'on a atteint T, vaut 0 si aux n'est pas dans T, un entier positif strict. sinon
        tb=estdans(aux(1),aux(2),z); // on teste si l'on cree une boucle
        if (tT>0) then
            z=[z ; aux]; // si l'on a atteint T, on s'arrete    
        elseif (tb>0) then 
            z=z(1:1:tb,:) // on efface la boucle  
        elseif (tb==0) then
            z=[z ; aux];
        end
    end
    res=z;
endfunction

// genere et trace une excursion de marche a boucles effacees issue de (a,b) dans la grille n x n jusqu'a ce qu'elle touche T
function res=LERWplot(a,b,T,n)
    z=[a b];
    tb=0;
    tT=0;
    isoview(0,n-1,0,n-1);
    while (tT==0)
        aux=whatsnext(z($,1),z($,2),n); // prochain point de la marche
        tT=estdans(aux(1),aux(2),T); // on teste si l'on a atteint T, vaut 0 si aux n'est pas dans T, un entier positif strict. sinon
        tb=estdans(aux(1),aux(2),z); // on teste si l'on cree une boucle
        if (tT>0) then
            isoview(0,n-1,0,n-1);
            xset('thickness',2) 
            plot2d([z($,1) aux(1)],[z($,2) aux(2)],style=5) // on trace le dernier pas
            z=[z ; aux]; // si l'on a atteint T, on s'arrete    
        elseif (tb>0) then 
            z=z(1:1:tb,:) // on efface la boucle  
            clf();
            isoview(0,n-1,0,n-1);
            xset('thickness',2) 
            plot2d(z(:,1),z(:,2), style=5) // on trace le chemin sans la boucle
        elseif (tb==0) then
            isoview(0,n-1,0,n-1);
            xset('thickness',2) 
            plot2d([z($,1) aux(1)],[z($,2) aux(2)],style=5) // on trace le dernier pas
            z=[z ; aux];
        end
    end
    res=z;
endfunction


// la fonction qui genere l'arbre uniforme
function UST(n)
    T=[2 3];
    a=0;
    for i=0:1:n-1
        for j=0:1:n-1
            test=estdans(i,j,T);
            if (test==0) then
                 branche=LERW(i,j,T,n);
                 xset('thickness',2)
                 plot2d(branche(:,1),branche(:,2),style=5)
                 T=[T ; branche];  
            else 
                a=a+1;      
            end   
        end
    end   
endfunction