/////////////////////////////////////////////////////////////////////////
//  Problème du voyageur de commerce 
//  Recuit simulé "optimisé" 
//  Explorations/propositions de duree d=nb de villes => tous les schéma de décroissance 
//  de température conviennent.
//  Lorsque d=1, l'algo est équivalent au recuit classique 
//
//
/////////////////////////////////////////////////////////////////////////
//Quelques valeurs : (T=100 ou 1000 ou 2000 ou 300) (Temp0=1 a 20 suivant schéma)
//(Nville=5, 10, 15 ou 20)
/////////////////////////////////////////////////////////////////////////
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 '); 

//T=5000;
a=.9;
//Nville=10;
//Temp0=1.5*sqrt(2)*Nville;
//Temp0=200;


//Placement aléatoire des villes sur [0,1]^2
m=Nville;
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 du circuit initial  
Sigma0=grand(1,'prm',[1:Nville]');
X0=Villes(Sigma0,:);
X0=[X0;Villes(Sigma0(1),:)];

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



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




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


for t=1:T

//Proposition d'une inversion "Sigmarev" d'une séquence de villes entre I et J

            Sigmarev=Sigma(:,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 de l'inversion proposée
if Urev<= U(t) then X=[X,Xrev]; 
                    U=[U,Urev];
                    Sigma=[Sigma,Sigmarev];
               else  if epsi(t)<= exp(-(Temp(t))^(-1)*(Urev-U(t))) then X=[X,Xrev]; 
                          U=[U,Urev];
                          Sigma=[Sigma,Sigmarev];
                          else X=[X,X(:,2*t-1:2*t)]; 
                               U=[U,U(t)];
                               Sigma=[Sigma,Sigma(:,t)];
                               end
              
 
               end
               
               

//Mise en mémoire du plus court circuit exploré 
if  U(t+1)<= Uopt then
Xopt=[X(:,2*t+1),X(:,2*t+2)];
Uopt=U(t+1);
Topt=t+1;
end               
 
 
//Schéma de décroissance de la température
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)];

//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+1),X(:,2*t+2),2)
//xtitle("Circuits selectionnes au cours du temps, nb de villes : "+string(Nville),"","")





end



xset("window",3);
xbasc();
xselect();
plot2d(1:T+1,U)
xtitle(" Variations des longueurs des circuits, nb villes : "+string(Nville),"axe du temps","longueur")
 
xset("window",4);
xbasc();
xselect();
plot2d(1:T+1,Temp,rect=[0,0,T+1,1])
xtitle("Decroissance de la temperature : ","axe du temps","")

xset("window",5);
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+1),X(:,(2*T+2)),2)
xtitle("Circuit terminal, nb d iterations : "+string(T),"","")



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),5)
xtitle("Circuit optimal visite par l algorithme : "+string(Topt),"","")