//simulations en amphi : a) cercle à 20 points, b) aleatoire unif c) étoile 8 branches d) étoile 25 branches

clear;
stacksize(150000000);
driver("X11");
clf(1);
scf(1);
//xset("wpos",600,0)
xset("wpos",400,0)
xselect();
f=gcf();
f.pixmap="on";
//xset("wdim",1000,1000);
xset("wdim",800,800);

clf(2);
scf(2);
//xset("wpos",600,0)
xset("wpos",400,0)
xselect();
f2=gcf();
f2.pixmap="on";
//xset("wdim",1000,1000);
xset("wdim",800,800);



Nvilles=35;
Niterations=0;
couleur=2;
scf(1);



places_villes=x_choose(["Aleatoire gaussien"; "Aleatoire uniforme";"Reseau regulier";"Etoile";"Cercle";"Aleatoire circulaire";"Exemple simple"],"Type de placement des villes","Quitter");

xselect();

select places_villes

	case 0 then break;
	
	case 1 then 
		
		rand("normal"); Villes=rand(Nvilles,2);  rand("uniform");
		absminaux=min(Villes(:,1));
		absmaxaux=max(Villes(:,1));
		ordminaux=min(Villes(:,2));
		ordmaxaux=max(Villes(:,2));
		rect_data_bounds=1.005*[absminaux,ordminaux;absmaxaux,ordmaxaux];
		
	case 2 then 
		
		rand("uniform"); Villes=rand(Nvilles,2); Rmax=1.1*max(abs(Villes));rect_data_bounds=[0,0;1,1]; 
		
		
	case 3 then
		
		vecteur_quadrillage=[3,4,5,7,10,12,15,20,25,40,60];
		string_quad=[];
		for i=1:length(vecteur_quadrillage)
			string_quad=[string_quad;string(vecteur_quadrillage(i))];
		end
		
		rNvilles=vecteur_quadrillage(x_choose(string_quad,"Nombre de villes par ligne","Quitter"));
			xselect();
		Nvilles=rNvilles^2;
		Villes_abs=1:rNvilles;
		for j=2:rNvilles
			Villes_abs=[Villes_abs,1:rNvilles];
		end
		
		Villes_ord=ones(1,rNvilles);
		for j=2:rNvilles
			Villes_ord=[Villes_ord,j*ones(1,rNvilles)];
		end
		
		perm=grand(1,'prm',[1:Nvilles]');
		Villes=zeros(Nvilles,2);
		
		for i=1:Nvilles
			Villes(i,:)=[Villes_abs(perm(i)),Villes_ord(perm(i))];
		end
		
		Villes=Villes./rNvilles;
		
		rect_data_bounds=[0,0;1+1/rNvilles,1+1/rNvilles];
		
	case 4 then

		vecteur_branches=[5,6,8,10,12,15,20,25,40,60];
		string_quad=[];
		for i=1:length(vecteur_branches)
			string_quad=[string_quad;string(vecteur_branches(i))];
		end
		
		Nbranches=vecteur_branches(x_choose(string_quad,"Nombre de branches","Quitter"));
		xselect();
		Nvilles=2*Nbranches;
		Rmin=2;
		Rmax=3+2*floor(Nbranches/15);
		thetaint=linspace(0,2*%pi,Nbranches+1);
		thetaint=thetaint(1:$-1);
		thetaext=thetaint+%pi/Nbranches;
//		Villes_abs=[Rmin*cos(thetaint),Rmax*cos(thetaext)];
//		Villes_ord=[Rmin*sin(thetaint),Rmax*sin(thetaext)];
		
		Villes_etoile=[];
		for i=1:Nbranches
			Villes_etoile=[Villes_etoile;...
			Rmin*cos(thetaint(i)),Rmin*sin(thetaint(i));...
			Rmax*cos(thetaext(i)),Rmax*sin(thetaext(i))];
		end
		
//		rect_data_bounds=1.01*[-Rmax,-Rmax;Rmax,Rmax];
//		xclick();
//		xpoly(Villes_etoile(:,1),Villes_etoile(:,2));
//		q=get("hdl");
//		q.closed="on";
//		q.thickness=2;
//		q.foreground=couleur;
//		xselect();
//		show_pixmap();
//		xselect();
//		xclick();
//		q.data=[];
		
		perm=grand(1,'prm',[1:Nvilles]');
		Villes=[];
		
		for i=1:Nvilles
			Villes=[Villes;Villes_etoile(perm(i),:)];
		end
		
		rect_data_bounds=1.01*[-Rmax,-Rmax;Rmax,Rmax];
		
	case 5 then
	
		vecteur_villes=[5,6,8,10,12,15,20,25,40,60];
		string_quad=[];
		for i=1:length(vecteur_villes)
			string_quad=[string_quad;string(vecteur_villes(i))];
		end
		
		Nvilles=vecteur_villes(x_choose(string_quad,"Nombre de villes","Quitter"));
	xselect();
	
		theta=linspace(0,2*%pi,Nvilles+1);
		theta=theta(1:$-1);	
		Villes_abs=cos(theta);
		Villes_ord=sin(theta);
		
		perm=grand(1,'prm',[1:Nvilles]');
		Villes=zeros(Nvilles,2);
		
		for i=1:Nvilles
			Villes(i,:)=[Villes_abs(perm(i)),Villes_ord(perm(i))];
		end
		
		rect_data_bounds=[-1.2,-1.2;1.2,1.2];
		
	case 6 then 	
	
		rand("normal");
		
		theta=linspace(0,2*%pi,Nvilles+1);
		theta=theta(1:$-1)+rand(1,Nvilles)./40;
		Villes_abs=cos(theta).*(1+rand(1,Nvilles)/4);
		Villes_ord=sin(theta).*(1+rand(1,Nvilles)/4);
		
		rand("uniform");
		perm=grand(1,'prm',[1:Nvilles]');
		Villes=zeros(Nvilles,2);
		
		for i=1:Nvilles
			Villes(i,:)=[Villes_abs(perm(i)),Villes_ord(perm(i))];
		end
		
		absminaux=min(Villes(:,1));
		absmaxaux=max(Villes(:,1));
		ordminaux=min(Villes(:,2));
		ordmaxaux=max(Villes(:,2));
		rect_data_bounds=1.005*[absminaux,ordminaux;absmaxaux,ordmaxaux];

	
	else
		h=1;
		Nvilles=6;
		Villes=[0,0;0,1;-2,2.5;-1,1;2,-1.5;3,2;-3,-0.5];
		absminaux=min(Villes(:,1));
		absmaxaux=max(Villes(:,1));
		ordminaux=min(Villes(:,2));
		ordmaxaux=max(Villes(:,2));
		rect_data_bounds=1.005*[absminaux,ordminaux;absmaxaux,ordmaxaux];
		
end //pour select places_villes



h=300*(1+5*floor(Nvilles/50));
deff('y=temperature(t)','y=h*ones(t)./(t+1)');



Villesvoyage=Villes;



//for i=1:Nvilles
//	xstring(Villes(i,1)+.1,Villes(i,2)+.1,string(i));
//end


distances=zeros(Nvilles,Nvilles);
for i=1:Nvilles
	for j=1:Nvilles
	distances(i,j)=sqrt((Villes(i,1)-Villes(j,1))^2+(Villes(i,2)-Villes(j,2))^2);	
	end
end

function e=energie(ordre_aux)
	e=distances(ordre_aux(Nvilles),ordre_aux(1));
	for i=1:Nvilles-1
		e=e+distances(ordre_aux(i), ordre_aux(i+1));
	end
endfunction


function [i,j]=deuxvilles()
	x=grand(1,2,'uin',1,Nvilles);
	while x(1)==x(2)
		x=grand(1,2,'uin',1,Nvilles);
	end
	i=min(x);
	j=max(x);
endfunction







Rmax=1.1*max(abs(Villes));
xpoly(Villes(:,1),Villes(:,2));
p=get("hdl");
p.line_mode="off";
p.thickness=2;
p.mark_style=9;
p.mark_foreground=5;
p.mark_background=4;
p.mark_size=2;

a=gca();
a.data_bounds=rect_data_bounds;
a.title.text="Cliquer pour afficher le parcours initial";
a.title.font_size = 5; 
a.title.font_style = 6; 
show_pixmap();
xselect();

 xclick();
a.title.text="Le voyageur de commerce par le recuit-simule";
xpoly(Villesvoyage(:,1),Villesvoyage(:,2));
q=get("hdl");
q.closed="on";
q.thickness=2;
q.foreground=couleur;
xselect();
show_pixmap();

ordre=1:Nvilles;
energie_actuelle=energie(ordre);
suite_energie=[energie_actuelle];

vect_choix=[20,500,1000,7000,10000,15000];

string_choix=[];
for i=1:length(vect_choix)
	string_choix=[string_choix;"+"+string(vect_choix(i))];
end

while %T
	Choix=x_choose([string_choix;"Longueurs des parcours"],string(Niterations)+" iterations de l''algorithme","Quitter");
	
	
	
	select Choix
	
	case 0 then  break;
	
	case length(vect_choix)+1 then
		
		clf(2);
		scf(2);
		abs_diag=0:Niterations;
		plot2d(0,0)
		plot2d(abs_diag,suite_energie,style=5,axesflag=5)
		b=gca();
		b.title.text="Longueur du parcours en fonction du nombre d''iterations";
		b.title.font_size = 5; 
		b.title.font_style = 6; 
		b.x_label.font_size = 5; 
		b.x_label.font_style = 6; 
		b.x_label.text="Nombre d''iterations de l''algorithme";
		b.y_label.font_size = 5; 
		b.y_label.font_style =6; 
		b.y_label.text="Longueur du parcours du voyageur";				
		show_pixmap();
		xselect();
		
		clf(1);
		scf(1);
		xpoly(Villes(:,1),Villes(:,2));
		p=get("hdl");
		p.line_mode="off";
		p.thickness=2;
		p.mark_style=9;
		p.mark_foreground=5;
		p.mark_background=4;
		p.mark_size=2;
		
		a=gca();
		a.data_bounds=rect_data_bounds;
		a.title.text="Le voyageur de commerce par le recuit-simule";
		a.title.font_size = 5; 
		a.title.font_style = 6;
		xpoly(Villesvoyage(:,1),Villesvoyage(:,2));
		q=get("hdl");
		q.closed="on";
		q.thickness=2;
		q.foreground=couleur;
		
		
	else 
	
		for t=1:vect_choix(Choix)
			Temp=temperature(Niterations);
			[x,y]=deuxvilles();
			if x>1 then
				if y<Nvilles then
					ordre_candidat=[ordre(1:x-1),ordre(y:-1:x),ordre(y+1:Nvilles)];
				else
					ordre_candidat=[ordre(1:x-1),ordre(Nvilles:-1:x)];
				end //pour 	if x(2)<Nvilles then
			else
				if y<Nvilles then
					ordre_candidat=[ordre(y:-1:1),ordre(y+1:Nvilles)];
				else
					ordre_candidat=[ordre(Nvilles:-1:1)];
				end //pour 	if x(2)<Nvilles then		
			end //pour if x(1)>1 then
			energie_candidat=energie(ordre_candidat);
			
			if (energie_candidat<energie_actuelle|rand(1,1)<exp((energie_actuelle-energie_candidat)/Temp)) then
				ordre=ordre_candidat;
				energie_actuelle=energie_candidat;
				suite_energie=[suite_energie,energie_actuelle];
				for v=1:Nvilles
					Villesvoyage(v,:)=Villes(ordre(v),:);
				end // pour for v=1:Nvilles
				q.data=Villesvoyage;
				xselect();
				show_pixmap();

			else
				suite_energie=[suite_energie,energie_actuelle];
				//affrate=ordre
			end // pour if (energie_candidat<energie_actuelle|rand(1,1)<exp((energie_actuelle-energie_candidat)/Temperature)) then
			Niterations=Niterations+1;		
		end //pour for t=1:vect_choix(Choix)
	end //pour select Choix
end //pour while %T

f.pixmap="off";
f2.pixmap="off";
driver("Rec");