function[pi]=tr_env_auc_cor_unif(cardinal_alea)
// no correlation + uniformity
// environment = i.i.d. with uniform common law
  dd=cardinal_alea
  pi=1/dd*ones(dd,dd)
endfunction

function[pi]=tr_env_cor_proche_vois(cardinal_alea)
// correlation with closest neighbours
  dd=cardinal_alea
  pi=diag([0.5 1/3*ones(1:(dd-2)) 0.5])+ diag([0.5 1/3*ones(1:(dd-2))],1)+...
       diag([1/3*ones(1:(dd-2)) 0.5],-1)
endfunction

function[pi]=tr_env_cor_crois(cardinal_alea,rho)
// correlation with increasing trend
// rho is the correlation intensity 0.8
  dd=cardinal_alea
  pi=diag([rho*ones(1:(dd-1)) 1])+diag([(1-rho)*ones(1:(dd-1))],1)
endfunction

function[pi]=tr_env_cor_decrois(cardinal_alea,rho)
// correlation with decreasing trend
// rho is the correlation intensity 0.8
  dd=cardinal_alea
  pi=diag([1 rho*ones(1:(dd-1))])+diag([(1-rho)*ones(1:(dd-1))],-1)
endfunction

function[pi]=tr_env_cor_et_chute(cardinal_alea,rho)
// correlation with probability of falling
// rho is the correlation intensity 0.9
  dd=cardinal_alea
  pi=diag([rho*ones(1:dd)])
  pi(2:$,1)=1-rho
  pi(1,2)=1-rho
endfunction

function[pi]=tr_env_cor_et_ascension(cardinal_alea,rho)
// correlation with a probability of rise
// rho is the correlation intensity 0.9
  dd=cardinal_alea
  pi=diag([rho*ones(1:dd)])
  pi(1:$,$)=1-rho
  pi(dd,dd-1)=1-rho
endfunction



function[Hypermatrices]=constr_HyperM(taille,env,control,pi,...
controlled_dynamics) 
// Construction of a family of matrices of transition probabilities on 
//  \{1,...,n\} X \{1,...,m\} from the problem data:
// Hypermatrices is a list of hypermatrices indexed by the control
// control is a vector
// controlled_dynamics(x,r,u) is a function with real values
dims=[prod(size(taille)),prod(size(env))];
ddims=[dims(1),dims(1),dims(2),dims(2)];
Hypermatrices=list();

cardinal_control=prod(size(control));
cardinal_taille=prod(size(taille));
for l=1:cardinal_control

//Hypermat=hypermat([ddims],sparse([],[],[prod(dims),prod(dims)]));
Hypermat=hypermat([ddims]);
// allows a natural expression for transitions from a product space
// towards itself


for y=1:dims(2)
image=controlled_dynamics(taille,y,control(l));
// vector of the images of the vector state for r fixed

[indexes_image_discretisee]=predec_succes(taille,image);
indexes1=indexes_image_discretisee(1);
indexes2=indexes_image_discretisee(2);
probabilites=indexes_image_discretisee(3);

M1=zeros(cardinal_taille,cardinal_taille);
M2=zeros(M1);
for i=1:cardinal_taille,
  M1(i,indexes1(i))=probabilites(i);
  M1(i,indexes2(i))=1-probabilites(i);
end

for z=1:dims(2)
Hypermat(:,1,y,z)=(1-ateb)*pi(y,z);
// Alive (index>1), the plant faces a death probability
// and the environment transitions are those of the matrix pi 
 Hypermat(:,:,y,z)=Hypermat(:,:,y,z) + (M1+M2)*ateb*pi(y,z);
end

end
// and not Hypermat(x,dyn,y,:)=(1-ateb)*pi(y,:) end, end
// because, at the previous step, one may have attributed
// Hypermat(x,dyn,y,:) with Hypermat(x,1,y,:)=ateb*pi.
Hypermatrices(l)=Hypermat
end
endfunction



function[Matrices]=conv_HyperM(Hypermatrices)
//Matrices is a list of transition matrices indexed by
//the control
dd=size(Hypermatrices(1))
dims=[dd(1),dd(3)]
//dimensions of the hypermatrix
Matrices=list()
for i=1:prod(size(control))
Mat=[]
for j=1:dims(2)
if dims(2)<>1
 a=Hypermatrices(i)(:,:,j,:).entries
 b=matrix(a,dims(1),prod(dims))
else
 b=matrix(Hypermatrices(i)(:,:,j,:),dims(1),prod(dims))
end
Mat=[Mat;b]
// Transformation into transition matrix (prod(dims),prod(dims))
// This latter matrix M is adapted to the resolution of the SDPE
end
Matrices(i)=full(Mat)
end
endfunction


function [n]=convert1(i,j,dims)
// i and j are coordinates in a hypermatrix of dimension dims
// n is the position of (i,j) in a vector of dimension prod(dims)
if i>dims(1) | j>dims(2)
 n='i ou j > dims'
else
 n=(j-1)*dims(1)+i;
end
endfunction


function [i,j]=convert2(n,dims)
// n is the position of (i,j) in a vector of dimension prod(dims)
// i and j are scoordinates in a hypermatrix of dimension dims
    i=modulo(n,dims(1))
    ind=find(i==0)
    i(ind)=dims(1)
    j=int(n/dims(1))+1
    ind2=find((dims(1))\n==int((dims(1))\n))
    jj=(dims(1))\n
    j(ind2)=jj(ind2)
endfunction


function[Hmat]=convert2_mat(matrix,dims)
// matrix is a matrix (ex:feed) whose elements come from 
// a vector space with dimension greater than 2
// Hmat is the hypermatrix corresponding to this matrix (ex:feed_hyp)
  Hmat=hypermat([dims,T-1])
  n=1:prod(dims)
  [i,j]=convert2(n,dims)
  for t=1:T-1
  for index=1:prod(dims)
  Hmat(i(index),j(index),t)=matrix(n(index),t)
  end
  end
endfunction


function[instant_cost]=conv_cost(taille,env,control,cost,horizon)
//instant_cost is a list of instantaneous costs indexed for 
// each value of the control
//taille is a vector
//env is a vector
//control is a vector
//cost is a function(x,u,t)
dims=[prod(size(taille)),prod(size(env))]
instant_cost=list()
for j=1:prod(size(control))
cost_i=[]  
  for t=1:horizon
// cost_it = hypermat(dims,sparse([],[],[dims(1),dims(2)]));
cost_it = hypermat(dims);
for state1=1:dims(1)
for state2=1:dims(2)
  cost_it(state1,state2)=cost(taille(state1),env(state2),control(j),t)
end 
end
cost_it=matrix(cost_it,1,prod(dims))
// Transformation into vector : size(cost_i)=[1,prod(dims)]
// adapted to the resolution of the SDPE
cost_i=[cost_i (cost_it)']
  end
instant_cost(j)=full(cost_i)
end
endfunction