// exec('plantII.sce');

getf('plantI2.sci');
getf('plantII.sci');

T=10;
ateb=0.9;
rr=1.2;
power = 0.5;
xp=(ateb*rr*power)^{1/(1-power)};
xm=(xp/rr)^{1/power};
taille=[0:(xm/3):(1.2*xp)];
//taille=0:2;
env=[0.6 0.8 1 1.2 1.4];
// env=[10:12]
control=0:0.05:1;
// control=0:0.5:1;

cardinal_env=prod(size(env));
cardinal_taille=prod(size(taille));

[pi]=tr_env_auc_cor_unif(cardinal_env);  // no correlation + uniformity
//[pi]=tr_env_cor_proche_vois(cardinal_env); //cor with close neighbours
//[pi]=tr_env_cor_crois(cardinal_env,rho); 
//cor (rho)  with trend to grow
//[pi]=tr_env_cor_decrois(cardinal_env,rho); 
//cor (rho) with trend to decrease
//[pi]=tr_env_cor_et_chute(cardinal_env,rho); 
//cor (rho) with a probability of fall
//[pi]=tr_env_cor_et_ascension(cardinal_env,rho); 
//cor (rho) with a probability of rise
//[pi]=tr_env_cor_crois(cardinal_env,1); //constant env: total correlation

function[b]=dyn_plant_control(x,r,v)
//b = total biomasse
//x = vegetative biomass
//r = environment
//v = control
     b=v*r*x^{power}
endfunction

controlled_dynamics=dyn_plant_control;

[Hypermatrices]=constr_HyperM(taille,env,control,pi,controlled_dynamics);

[transition_matrix]=conv_HyperM(Hypermatrices);


function[ci]=offspring2(taille,env,control,time)
// fitness of the plant
  ci=(1-control)'*env*taille^(power)*ateb^{time-1}
  ind=find(taille==1),ci(ind)=0
endfunction

[offspring]=conv_cost(taille,env,control,offspring2,T);

final_cost_nul=zeros( cardinal_env * cardinal_taille ,1); 
// the final cost is zero



dims=[cardinal_taille,cardinal_env];

stacksize(3000000);

instant_cost=offspring;
final_cost=final_cost_nul;

[value,feedback]=Bell_stoch(transition_matrix,instant_cost,final_cost,1);

index_taille_init=grand(1,1,'uin',2,cardinal_taille);
index_env_init=grand(1,1,'uin',1,cardinal_env);

[initial_state]=convert1(index_taille_init,index_env_init,dims);

[z]=trajopt(transition_matrix,feedback,instant_cost,final_cost,initial_state);
    
// computation of optimal trajectories (indexes)

[index_taille,index_env]=convert2(z(1),dims)

x=taille(index_taille);
e=env(index_env);
// optimal state trajectories in naturel units

xset("window",1);xbasc();plot2d2(1:prod(size(x)),x);
xtitle("taille")

xset("window",2);xbasc();plot2d2(1:prod(size(e)),e,rect=[0,0,T+1,max(e)+1]);
xtitle("environment")

xset("window",3);xbasc();plot2d2(1:prod(size(control(z(2)))),control(z(2)));
xtitle("control")

xset("window",4);xbasc();plot2d2(1:prod(size(z(3))),z(3));
xtitle("fitness")