// \begin{verbatim}
// exec proba_extinction.sce

// demographic stochasticity parameters
sig_e=0.2; sig_d=0.1; rbar=0.03;

// environmental stochasticity parameters
w_min=-1;w_max=1;p=0.5;q=1-p;

// Ricker dynamics parameters 
Cap=0; 			
// carrying capacity

// constraints thresholds
N_min=1; N_max=10^3;
h_min=0.2;h_max=0.2;

// time horizon
Horizon=250 ;


function y=dynpop(x,w)
  y=zeros(x);
  z=x(x>=N_min);
  y(x>=N_min)= z.*(1+rbar+((sig_e^2+sig_d^2 ./z)^0.5).* w(x>=N_min));
// zero when x < N_min
endfunction

function y=dyna_exp(x,h,w)  
	y=dynpop(x-h,w); 				
endfunction

// Indicator function of state constraint
function y=Ind_x(t,x) 
  y=bool2s(x>=N_min);
endfunction

// Characteristic function of control constraint (approximate)
function [y]=Phi_u(t,x,h) 
  SGN = bool2s(h_min <= h) ; 
  y=0*SGN + 1/%eps *(1-SGN) ;
endfunction



/////////////////////////////////////
// State and control Discretization
/////////////////////////////////////

// Grid state x
x_min=0; x_max=N_max;delta_x=1;
grille_x=x_min:delta_x:x_max;
S=size(grille_x);
NN=S(2);

// Grid control 
u_min=0; u_max=h_max;delta_u=0.05;
u=u_min:delta_u:u_max;
R=size(u); MM=R(2);

function [z]=Projbis(x)
	z=round(x./delta_x).*delta_x;
	z=min(z,x_max);
	z=max(z,x_min);
endfunction

function i=Indice(x)
	i=int((x-x_min)./delta_x)+1;
endfunction

// Discretized dynamics
function [z]=dyna(x,h,w)
	xsuiv=dynpop(x-h,w);
	P=Projbis(xsuiv);
	z=Indice(P);
	endfunction 

for (ii=1:NN)
	Etat_x(ii)=x_min+(ii-1)*delta_x;
end;

for (ii=1:MM)
	Control_u(ii)=u_min+(ii-1)*delta_u;
end;

//////////////////
// Graphics
///////////////////

xset("window",1);xbasc(1);
xtitle("Maximal viability probability",...
	   'x','Max_{h}P(N(T)>=N_min)');

xset("window",2); xbasc(2);
xtitle("Population",'time t','N(t)');

xset("window",3); xbasc(3);
xtitle("Catches",'time t','h(t)');


///////////////////////	
//	  Dynamic programming
///////////////////////

x=grille_x;
W=zeros(Horizon,NN);

// Initialization at horizon T

for (i=1:NN) 
	W(Horizon,i)=Ind_x(Horizon,Etat_x(i)); 
	end

// Bellman equation
for(t=Horizon-1:-1:1)
	for (i=1:NN)
		xx=Etat_x(i);
		for (j=MM:-1:1)
			uu=Control_u(j);
			g_min(j)=-Phi_u(t,xx,uu)+W(t+1,dyna(xx,uu,w_min)); 
			g_max(j)=-Phi_u(t,xx,uu)+W(t+1,dyna(xx,uu,w_max)); 
		  g(j)=p*g_min(j)+(1-p)*g_max(j); // expected value
		end,
		[Vopt,jopt]=max(g) ; 
		W(t,i)=Vopt;
		j_opt(t,i)=jopt;
	end,
end,



// Viability kernel
Viab=W(1,:)';K=1-W(Horizon,:)';
xset("window",1);xbasc(1);
plot2d(Etat_x,Viab,-4,rect=[0,0,x_max,1]);
xtitle("Maximal probability viability","initial abundance",...
	   "probability");
// legends("Maximal probability viability",-4,"lr");
 
 
////////////////////////////////////////////////////
// Simulations 
////////////////////////////////////////////////////

// N0=200;
threshold_viab=0.8

if (max(W(1,:)) >= threshold_viab) then 
for (k=1:10)
	i(1)=int(rand(1)*(NN-1))+1;compt=0;
	while(W(1,i(1))<threshold_viab ) 
	then i(1)=int(rand(1)*(NN-1))+1;end, 
	// find viable state x(i(1))
	traj_x(1)=Etat_x(i(1));
	//traj_x(1)=N0; i(1)=Indice(N0);
	for (t=1:Horizon-1)
		if (rand(1)<p) then w(t)=w_min; else w(t)=w_max; end,
		h_viab(t)=Control_u(j_opt(t,i(t)));
		i(t+1)=dyna(traj_x(t),h_viab(t),w(t));
		traj_x(t+1)=Etat_x(i(t+1));
	end,
	Tempx=(1:Horizon)-1;
	Nmin=zeros(1,Horizon)'+N_min;
	Nmax=zeros(1,Horizon)'+N_max;
	NLin=traj_x(1)*(1+rbar).^Tempx';
	//
	xset("window",2);plot2d(Tempx,[traj_x Nmin Nmax NLin],...
	[1,2,3,4],leg="N(t)@Nmin@Nmax@Ncertain",...
							rect=[1,N_min,Horizon,N_max]);
	xset("window",2);plot2d(Tempx,[traj_x Nmin Nmax NLin],...
							rect=[1,N_min,Horizon,N_max]);
	plot2d(Tempx,[traj_x Nmin Nmax NLin],style=-[1,2,3,4])
	legends(["N(t)";"Nmin";"Nmax";"Ncertain"],-[1,2,3,4])
	// STOP...
	Temps=1:Horizon-1;
	xset("window",3);plot2d(Temps,h_viab,...
							rect=[1,u_min,Horizon,u_max]);

end,
else " beta kernel is empty "	
end

// \end{verbatim}