//Conditions de bonne utilisation : commencer par réduire la console scilab de façon à ce que elle ne chevauche pas la fenetre du dessous (que l'on ne met jamais au premier plan, pour des raisons de manque de surface sur l'écran).

clear;
stacksize(150000000);
driver("X11");
rand("uniform");
//deff('y=changecouleur(t)','y=modulo(t+1,35)+1');

nbre_points_curve=1000;
esp_geometrique=4;
nbre_etapes_down=20;
coeff_duree=1; //pour une valeur 1, les parametres sont choisis de façon à ce que l'animation se déroule en un temps adapté à l'amphi
nbre_histo_fluct=floor(60*coeff_duree);
nbre_sim=0;

%symbol=1; 
%times=2; 
%times_italic=3; 
%times_bold=4; 
%helvetica= 6; 
%helvetica_italic=7; 
%helvetica_bold=8; 
%helvetica_bold_italique=9; 
%8pts=0;%10pts=1;%12pts=2; 
%14pts=3;%18pts=4;%24pts=5; 

a=1.7;
b=0.1;
deff('y=fct1(x)','y=(x).*(3-x).*(3-x).*((x-a).*(x-a)+b).*( ( ones(x)./( (x-2.6).*(x-2.6) +1 ) ) +ones(x) ).*( ( ones(x)./( (x-0.9).*(x-0.9) +1 ) ) +ones(x) )./(3.1-x)^2');//
C=1./intg(0,3,fct1);
deff('y=fct_densite(x)','y=C*fct1(x)');

echantillon=[];
ymax=esp_geometrique/3;
ymax2=1.05*fct_densite(2.74);
ymax3=1.001*fct_densite(2.74);
ymin=ymax/200;
rect_haut=[0,0,3,ymax];


x_factice=[0.5, 1, 2.5, 2, 0.3,1.7,2.3, 0.7, 2.95, 1.5,2.6];
nbre_factice=length(x_factice);
u=rand(1,nbre_factice);
u=(u-ones(u)./2)./100;
x_factice=x_factice+u;

u=rand(1,nbre_factice);
u=u./100;
y_factice=([2, 1.5, 1.1,1.1, 1, 6,1.2, 2,2.5,11,0.7]+u).*fct_densite(x_factice);
y_factice(5)=ymax*(0.9-rand(1,1)./10);

y_factice($-1)=min(y_factice($-1), 0.97*ymax);

abs_curve=linspace(0,3,nbre_points_curve);
ord_curve=fct_densite(abs_curve);

function style_haut() 
	a = gca(); // a est le handle sur ces axes par defaut 
	a.title.font_size = %24pts; //(=%18pts;)
	a.title.font_style = 6;  
	a.y_label.font_size = %18pts; //(=%14pts;)
	a.y_label.font_style =6;  
	a.title.text="Simulation par la methode du rejet";
	a.thickness=2; 
	xselect();
	show_pixmap();
endfunction 

function style_bas() 
	a = gca(); // a est le handle sur ces axes par defaut 
	a.title.font_size = %24pts; //(=%18pts;)
	a.title.font_style = 6; //(= %times_bold;)
	a.title.text="Histogramme des valeurs obtenues";
	a.thickness=2;
endfunction 


function premieraffichage()
	scf(1); 
	clf(1);
	style_haut();
	plot2d(abs_curve,ord_curve,style=5,rect=rect_haut,frameflag=1);
	legend("densite");
	xselect();
	show_pixmap();
endfunction 

function premieraffichageprime()
	scf(1); 
	clf(1);
	style_haut();
	plot2d(abs_curve,ord_curve,style=5,rect=rect_haut,frameflag=1);
	//legend("densite : y=x^2(3-x)((x-1.7)^2+0.1");
	xselect();
endfunction 
	
function premieraffichage_bas()
	scf(2); 
	clf(2);
	style_bas();
	plot2d(abs_curve,ord_curve,style=5,rect=rect_haut,frameflag=1);
	legend("densite");
	show_pixmap();
endfunction 

function affichage_bas2()
	scf(2); 
	clf(2);
	style_bas();
	plot2d(abs_curve,ord_curve,style=5,rect=rect_haut,frameflag=1);
	legend("densite");
endfunction 

//main()

scf(1); 
clf(1);
xset("wpos",180,0)
xselect();
f=gcf();
f.pixmap="on";
xset("wdim",820,300);

scf(2); 
clf(2);
xset("wpos",180,320)
//xselect();
f=gcf();
f.pixmap="on";
xset("wdim",820,330);

premieraffichage();

factice=%T;

Choix=1;
while (Choix<>0)
	Choix=x_choose(["Une simulation";...
				    "+1"; "+3"; "+5"; "+10";"+20";"+50";"+100 optimal";"+800 rapide";...
				    "Histo theorique";"Fluctus histo";... 
				    "Histo final"],string(nbre_sim)+" simulations","Quitter" );
	xselect();
	select Choix,
	
	case 1 then
		
		premieraffichageprime();
		if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
		show_pixmap();
		xclick();

		if factice then 
			for k=1:nbre_factice-1
				plot(x_factice(k),y_factice(k),"b+");
				xselect();
				show_pixmap(); 
				xclick();
			end
			x=x_factice($);
			y=y_factice($);
            plot(x,y,"b+");
			xselect();
			show_pixmap(); 
			pas_y=y/nbre_etapes_down;
			while y>=0
					 	yprime=y-pas_y;
					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
					 	y=yprime;
					 	show_pixmap();
			end
			plot(x,ymin,"r.");
			show_pixmap();	
			echantillon=[echantillon, x];
			nbre_sim=nbre_sim+1;	
			factice=%F;
		else
			bool=%T;
			while bool
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				if fct_densite(x)<y then 
					     plot(x,y,"b+");
						 xselect();
						 show_pixmap(); 
						 xclick();
					else plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
					 pas_y=y/nbre_etapes_down;
					 while y>=0
					 	yprime=y-pas_y;
					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
					 	y=yprime;
					 	show_pixmap();
					 end
					 plot(x,ymin,"r.");
					 show_pixmap();
					 bool=%F;
	
				end
			end
			echantillon=[echantillon, x];
			nbre_sim=nbre_sim+1;
		end
		
		case 2 then
		
		premieraffichageprime();
		if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
		show_pixmap();
		xclick();
		
		bool=%T;
		while bool
			x=3*rand(1,1);
			y=ymax*rand(1,1);
			if y>fct_densite(x) then 
					 plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
					 //xclick();
			else plot(x,y,"b+");
				 xselect();
				 show_pixmap(); 
				 pas_y=y/nbre_etapes_down;
				 while y>=0
				 	yprime=y-pas_y;
				 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
				 	y=yprime;
				 	show_pixmap();
				 end
				 plot(x,0,"r.");
				 show_pixmap();
				bool=%F;
			end
		end	
	echantillon=[echantillon, x];
	nbre_sim=nbre_sim+1;
	
	case 3 then
	
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
		
		for i=1:3
			bool=%T;
			while bool
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				if y>fct_densite(x) then 
						 plot(x,y,"b+");
						 xselect();
						 show_pixmap(); 
						 //xclick();
				else plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
//					 pas_y=y/nbre_etapes_down;
//					 while y>=0
//					 	yprime=y-pas_y;
//					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
//					 	y=yprime;
//					 	show_pixmap();
//					 end
					 plot(x,0,"r.");
					 show_pixmap();
					bool=%F;
				end
			end
		echantillon=[echantillon, x];
		end
		nbre_sim=nbre_sim+3;

	case 4 then
		
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
		
		for i=1:5
			bool=%T;
			while bool
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				if y>fct_densite(x) then 
						 plot(x,y,"b+");
						 xselect();
						 show_pixmap(); 
						 //xclick();
				else plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
//					 pas_y=y/nbre_etapes_down;
//					 while y>=0
//					 	yprime=y-pas_y;
//					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
//					 	y=yprime;
//					 	show_pixmap();
//					 end
					 plot(x,0,"r.");
					 show_pixmap();
					bool=%F;
				end
			end
		echantillon=[echantillon, x];
		end
		nbre_sim=nbre_sim+5;	

	case 5 then

	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();

		for i=1:10
			bool=%T;
			while bool
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				if y>fct_densite(x) then 
						 plot(x,y,"b+");
						 xselect();
						 show_pixmap(); 
						 //xclick();
				else plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
//					 pas_y=y/nbre_etapes_down;
//					 while y>=0
//					 	yprime=y-pas_y;
//					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
//					 	y=yprime;
//					 	show_pixmap();
//					 end
					 plot(x,0,"r.");
					 show_pixmap();
					bool=%F;
				end
			end
		echantillon=[echantillon, x];
		end
		nbre_sim=nbre_sim+10;
		
			
	case 6 then
	
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
		
		for i=1:20
			bool=%T;
			while bool
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				if y>fct_densite(x) then 
						 plot(x,y,"b+");
						 xselect();
						 show_pixmap(); 
						 //xclick();
				else plot(x,y,"b+");
					 xselect();
					 show_pixmap(); 
//					 pas_y=y/nbre_etapes_down;
//					 while y>=0
//					 	yprime=y-pas_y;
//					 	plot2d([x,x],[y, yprime],style=13,rect=rect_haut,frameflag=1);
//					 	y=yprime;
//					 	show_pixmap();
//					 end
					 plot(x,0,"r.");
					 show_pixmap();
					bool=%F;
				end
			end
		echantillon=[echantillon, x];
		end
		nbre_sim=nbre_sim+20;	
	
	case 7 then
		
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
		
		
		for i=1:50

			x=3*rand(1,1);
			y=ymax*rand(1,1);
			x_af=[x];
			y_af=[y];
			while  y>fct_densite(x) 
				x=3*rand(1,1);
				y=ymax*rand(1,1);
				x_af=[x_af,x];
				y_af=[y_af,y];
			end
			echantillon=[echantillon, x];
			plot(x_af,y_af,"b+");
			//plot2d([x,x],[y, 0],style=13,rect=rect_haut,frameflag=1);
			plot(x,0,"r.");
			xselect();
			show_pixmap();
		end
		nbre_sim=nbre_sim+50;	


	case 8 then
		
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
		
		abs2=linspace(0,3,50);
		plot(abs2,ymax2*ones(abs2),"y:");
		show_pixmap;
		xclick();
		
		for i=1:100
			
			x=3*rand(1,1);
			y=ymax2*rand(1,1);
			x_af=[x];
			y_af=[y];
			while  y>fct_densite(x) 
				x=3*rand(1,1);
				y=ymax2*rand(1,1);
				x_af=[x_af,x];
				y_af=[y_af,y];
			end
			echantillon=[echantillon, x];
			plot(x_af,y_af,"b+");
			//plot2d([x,x],[y, 0],style=13,rect=rect_haut,frameflag=1);
			plot(x,0,"r.");
			xselect();
			show_pixmap();
		end
		nbre_sim=nbre_sim+100;	
		
		case 9 then
		
	premieraffichageprime();
	if nbre_sim>0 then plot(echantillon,zeros(echantillon),"r.");  end
	show_pixmap();
	xclick();
	
		abs_rap=[];
		ord_rap=[];
		for i=1:800
			
			x=3*rand(1,1);
			y=ymax3*rand(1,1);
			while  y>fct_densite(x) 
				x=3*rand(1,1);
				y=ymax2*rand(1,1);
			end
			abs_rap=[abs_rap,x];
			ord_rap=[ord_rap,y];
			//plot2d([x,x],[y, 0],style=13,rect=rect_haut,frameflag=1);	
		end	
		nbre_sim=nbre_sim+800;
		echantillon=[echantillon, abs_rap];
		plot(abs_rap, ord_rap, "b+");
		plot(abs_rap,zeros(abs_rap),"r.");
		xselect();
		show_pixmap();
		
		case 10 then
		
			scf(2);
			clf();
			premieraffichage_bas();
			scf(1);
			
		case 11 then
			
			scf(2);
			clf();
			premieraffichage_bas();
			nbre_col_histo=min(nbre_sim,50);
			for j=1:nbre_histo_fluct
				affichage_bas2();
				n=max(floor(nbre_sim*j/nbre_histo_fluct),2);
				histplot(nbre_col_histo,echantillon(1:n),rect=rect_haut,frameflag=1);
				show_pixmap();
			end
			scf(1);
	
		case 12 then
			nbre_col_histo=min(nbre_sim,50);
			affichage_bas2();
			histplot(nbre_col_histo, echantillon,rect=rect_haut,frameflag=1);
			show_pixmap();
			scf(1);
	    
		
				
	else echantillon=[];	
		
	end
	xselect();
end	
scf(1);	
f=gcf();
f.pixmap="off";
scf(2);
f=gcf();
f.pixmap="off";
driver("Rec");