// exec('TP_plante_3.sce')
getf('TP_plante_2.sci');


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

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

[pi]=tr_env_auc_cor_unif(cardinal_env);  //aucune corrélation + uniformité
//[pi]=tr_env_cor_proche_vois(cardinal_env); //cor avec les proches voisins
//[pi]=tr_env_cor_crois(cardinal_env,rho); 
//cor (rho)  avec tendance à croitre
//[pi]=tr_env_cor_decrois(cardinal_env,rho); 
//cor (rho) avec tendance à décroitre
//[pi]=tr_env_cor_et_chute(cardinal_env,rho); 
//cor (rho) avec une probabilité de chute 
//[pi]=tr_env_cor_et_ascension(cardinal_env,rho); 
//cor (rho) avec une proba d'ascension
//[pi]=tr_env_cor_crois(cardinal_env,1); //env constant : cor totale

function[b]=dyn_plante_com(x,r,v)
//b = biomasse totale que la plante la plante peut allouer, sur
//    une saison, entre biomasse végétative et biomasse reproductive
//x = biomasse végétative
//v = la commande
     b=v*r*x^{puis}
endfunction

dynamique_commandee=dyn_plante_com;

[Hypermatrices]=constr_HyperM(taille,env,commande,pi,dynamique_commandee);

[matrice_transition]=conv_HyperM(Hypermatrices);


function[ci]=rejetons2(taille,env,commande,temps)
//Fonction définissant la fitness de la plante
  ci=(1-commande)'*env*taille^(puis)*beta^{temps-1}
  ind=find(taille==1),ci(ind)=0
endfunction

[rejetons]=conv_cout(taille,env,commande,rejetons2,T);

cout_fin_nul=zeros( cardinal_env * cardinal_taille ,1); 
//le cout final et nul

dims=[cardinal_taille,cardinal_env];

stacksize(3000000);

cout_instantane=rejetons;
cout_final=cout_fin_nul;

[valeur,feedback]=Bell_stoch(matrice_transition,cout_instantane,cout_final);

indice_taille_init=grand(1,1,'uin',2,cardinal_taille);
indice_env_init=grand(1,1,'uin',1,cardinal_env);

[etat_initial]=convert1(indice_taille_init,indice_env_init,dims);

[z]=trajopt(matrice_transition,feedback,cout_instantane,cout_final,etat_initial);
    
// calcul des trajectoires optimales (indices)

[indice_taille,indice_env]=convert2(z(1),dims)

x=taille(indice_taille);
e=env(indice_env);
// trajectoires optimales de l'état en unités naturelles

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("environnement")

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

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