/////////////////////////////////////////////////////////////////////////
//  Problème du voyageur de commerce 
//  Recuits simulés en interaction, Nombre de recuits=Nombre de particules
//  Les recuits interagissent à haute température, et deviennent progressivement
//  indépendants au cours du temps.
//  A population infinie, la loi d'une particule coincide avec la mesure de Gibbs cible.
//
//  Lorsque le nombre de particule =1, l'algo est équivalent au recuit classique 
//  Les explorations/propositions sont suffisamment longues pour que
//  tous les schémas de décroissance de température fonctionnent.
//
//
/////////////////////////////////////////////////////////////////////////
//Quelques valeurs : (T=100 ou 1000 ou 2000 ou 300) (Temp0=1 a 20 suivant schéma)
//(Nville=5, 10, 15 ou 20), (Npart= de 1 a 10)
/////////////////////////////////////////////////////////////////////////
clear

T=input('Horizon temporel ? \n ') ;  
Nville=input('Nombre de villes ? \n ') ;
Temp0=input('Temperature initiale  ? \n ') ;
Choix=input('Choix des decroissances de temperature : \n 1) Logarithmique \n 2) lineaire  \n 3) Polynomial \n 4) exponentiel  \n '); 
NPart=input('Nombre de particules ? \n ') ;
m=Nville;



//Placement aléatoire des villes sur [0,1]^2
a=1;
Xville=grand(1,Nville,'def');
Yville=grand(1,Nville,'def');
Villes=[];
for n=1:Nville
Villes=[Villes;Xville(n),Yville(n)];
end


//Choix aléatoire des NPart circuits initiaux  
Sigma0=grand(NPart,'prm',[1:Nville]');
for i=1:NPart
P0(:,:,i)=Villes(Sigma0(:,i),:);
X0(:,:,i)=[P0(:,:,i); Villes(Sigma0(1,i),:)];
end


//Longueurs des circuits initiaux
U0=zeros(NPart,1);
for i=1:NPart
for n=2:Nville+1
U0(i)=U0(i)+sqrt((X0(n,1,i)-X0(n-1,1,i))^2+(X0(n,2,i)-X0(n-1,2,i))^2);
end
end


//Longueur maximale des circuits initiaux
UMSort=sort(U0);
IMopt=find(U0==UMSort(1));
XMopt=X0(:,:,IMopt);

xset("window",1);
xbasc();
xselect();
xset("thickness",2)
for n=1:Nville
plot2d(Xville(n),Yville(n),-1,rect=[0,0,1,1])
end
for i=1:NPart
plot2d(XMopt(:,1),XMopt(:,2),IMopt)
end
xtitle("Nb de villes : "+string(Nville),"","")





//Initialisation des paramètres
Sigma(:,:,1)=Sigma0;
CNRej=[];
U=U0;
epsi=grand(1,T,'def');
X=X0;
Temp=[Temp0];

//Topt=1;
//Xopt=X0(:,:,NPart);;
//Uopt=U0(NPart);


//Schémas de décroissance de la température
for t=1:T

if Choix==1
Temp=[Temp,Temp0/log(t+1)];
elseif Choix==2
Temp=[Temp,Temp0/(t+1)];
elseif Choix==3
Temp=[Temp,Temp0/(t+1)^a];
else Temp=[Temp,Temp0*exp(-a*t)];
end              

//Temp=[Temp,Temp0/log(t+1)];
//Temp=[Temp,Temp0/(t+1)^a];
//Temp=[Temp,Temp0*exp(-a*t)];
end



for t=1:T
for i=1:NPart

            //Proposition de NPart inversions "Sigmarev" de séquences de 
            //villes entre I et J

            Sigmarev=Sigma(:,i,t);
            //Exploration des inversions pendant un temps m
            for s=1:m
            Sigmaopt=Sigmarev;
            I=1+floor(Nville*grand(1,1,'def'));
            J=1+floor(Nville*grand(1,1,'def'));
            K=[I,J];
            K=-sort(-K);
            
               for p=0:(K(2)-K(1))
               Sigmaopt(K(1)+p)=Sigmarev(K(2)-p);
               end
            Sigmarev=Sigmaopt;
            end

Xrev=Villes(Sigmarev,:);
Xrev=[Xrev;Villes(Sigmarev(1),:)];

Urev=0;
for n=2:Nville+1
Urev=Urev+sqrt((Xrev(n,1)-Xrev(n-1,1))^2+(Xrev(n,2)-Xrev(n-1,2))^2);
end


//Acceptation ou rejet des inversions proposées

if Urev<= U(i,t) then X(:,:,i)=Xrev; 
                    U(i,1:t+1)=[U(i,1:t),Urev];
                    Sigma(:,i,t+1)=Sigmarev;
                 else  if epsi(t)<= exp(-(Temp(t))^(-1)*(Urev-U(i,t))) then  X(:,:,i)=Xrev; 
                                                                         U(i,1:t+1)=[U(i,1:t),Urev];
                                                                         Sigma(:,i,t+1)=Sigmarev;
                          else X(:,:,i)=X(:,:,i); 
                               U(i,1:t+1)=[U(i,1:t),U(i,t)];
                               Sigma(:,i,t+1)=Sigma(:,i,t);
                               end
              
 
                 end




//xset("window",2);
//xbasc();
//xselect();
//xset("thickness",2)
//for n=1:Nville
//plot2d(Xville(n),Yville(n),-1,rect=[0,0,1,1])
//end
//plot2d(X(:,2*t),X(:,2*t+2),2)
//xtitle("Nb de villes : "+string(Nville),"","")



end


V=(Temp(t+1)-Temp(t))*U(:,t);
G=exp(-V');
Gsup=sort(G);
G=G/Gsup(1);
P=G/sum(G);
F=cumsum(P);


E=rand(1,NPart,'uniform');
Ac=find(E<=G);
NAc=length(Ac);
Rej=find(E>G);
NRej=length(Rej);
CNRej=[CNRej,NRej];

if NRej>0 then Expo=-log(grand(1,NRej+1,'def'));
               Texpo=cumsum(Expo);  
               Vexp=Texpo/Texpo(NRej+1);
               Xmulti=[];
               j=1;
               k=1;       
               while j<=NRej, if Vexp(j)<=F(k) then Xmulti=[Xmulti,X(:,:,k)]; 
                                                            j=j+1;
                                                else k=k+1;
                                end                

               end
               h=1; l=1;
         while l<=NPart, if E(l)>G(l) then X(:,:,l)=Xmulti(h);
                               h=h+1; 
                               l=l+1;
                          else l=l+1;
                          end 
                            
         end
end


end


xset("window",3);
xbasc();
xselect();
for i=1:NPart
plot2d(1:T+1,U(i,1:T+1)',i)
end
xtitle("Energie de chaque configuration, Nb de particules : "+string(NPart),"","")

xset("window",4);
xbasc();
xselect();
plot2d(1:T+1,Temp)
xtitle("Evolution de la temperature : ","","")

xset("window",5);
xbasc();
xselect();
plot2d(1:T,CNRej)
xtitle("Nb de particules rejetees : "+string(NPart),"","")

USort=-sort(-U(:,T+1));
Iopt=find(U(:,T+1)==USort(1));
Xopt=X(:,:,Iopt);

xset("window",6);
xbasc();
xselect();
xset("thickness",2)
for n=1:Nville
plot2d(Xville(n),Yville(n),-1,rect=[0,0,1,1])
end
plot2d(Xopt(:,1),Xopt(:,2),Iopt)
xtitle("Trajet optimal, Nb de particules : "+string(NPart),"","")