// parameters
power=1; //  gamma
ateb=0.9  // survival probability
r=1.2;
T=10; 
x0=1;
// ATTENTION. The control is mathematically indexed by t=0,...,T-1.
// However, it is indexed by 1:T=1:horizon in Scilab
// The state is indexed by 1:(T+1)=1:(1+horizon) in Scilab

getf('plantI1.sci');
strategy=list();
strategy(1)=strategy_D;
strategy(2)=strategy_R;
strategy(3)=strategy_G;

correspondance=list();
correspondance(1)=string("death");
correspondance(2)=string("random");
correspondance(3)=string("growth");

// In what follows, we adopt the terminology 
// of the arguments of function traj_feedback
initial_state=1;
dynamics=dyn_plant_control;
alea=ceil(ateb-rand(1,T));
// T independsent realizations of a binomial law B(ateb,1)
inst_cost=offspring;
final_cost=final_cost_zero;

for i=1:3 do
[b,v]=traj_feedback...
          (initial_state,dynamics,strategy(i),alea);
[f]=total_cost(b,v,inst_cost,final_cost);

 xset("window",3*(i-1));xbasc();plot2d2(1:T+1,f,rect=[0,0,T+2,max(f)+1]);
 xtitle("strategy "+correspondance(i)+" : fitness");
 
 xset("window",3*(i-1)+1);xbasc();plot2d2(1:T,v,rect=[0,0,T+1,max(v)+1]);
 xtitle("strategy "+correspondance(i)+" : control");
 
 xset("window",3*(i-1)+2);xbasc();plot2d2(1:T+1,b,rect=[0,0,T+2,max(b)+1]);
 xtitle("strategy "+correspondance(i)+" : state"); 
 
 printf("the total fitness for strategy "...
        +correspondance(i)+" is : %"+" f\n", f($))
 
 halt();
 xdel((3*(i-1)):(3*(i-1)+2));
// deletes the windows
end


// parameters
power = 0.5;
xp=(ateb*r*power)^{1/(1-power)};
xm=(xp/r)^{1/power};


strategy(4)=strategy_O;
correspondance(4)=string("optimal");

P=0.1:0.2:0.9;
for j=1:prod(size(P)) do
  power=P(j);
  xp=(ateb*r*power)^{1/(1-power)};
  xm=(xp/r)^{1/power};
x0=xm/10;
// initial biomass less than threshold xm
printf("gamma=%"+" f\n", P(j));
     for i=1:4 do
     [y,v]=traj_feedback(x0,dynamics,strategy(i),alea);
     [f]=total_cost(y,v,inst_cost,final_cost);
     printf("the total fitness for strategy "...
        +correspondance(i)+" is : %"+" f\n", f($))
     end
end