//\begin{verbatim}
stacksize(2*10^8);

// -----------------------------------------
// LABELLING STATES AND DYNAMICS
// -----------------------------------------

MM=3^{NN} ; // number of states

Integers=(1:MM)' ;  // labelled states

State=zeros(MM,NN); 
// Will contain, for each integer z in Integers,
// a sequence s=(s_1,...,s_NN)  with 
// s_1 \in \{1,2\} 
// and s_k \in \{0,1,2\} for k \geq 2, such that 
// 1) z = s_1 + s_2*3^{1} + ... + s_NN*3^{NN-1}
// 2) a mine profile p_1,...,p_NN is given by 
// p_k = y_1 + ... + y_k where y_1= s_1-1 \in \{0,1\}
// and y_j = s_j - 1 \in \{-1,0,1\} for j>1.

Increments=zeros(MM,NN); 
// Will contain, for each integer z in Integers,
// the sequence y=(y_1,...,y_NN).

// The initial profile is supposed to be p(0)=(0,0,...,0) 
// to which corresponds y(0)=(0,0,...,0) and
// s(0)=(1,1,...,1) and z(0)=1+ 3^{1} + ... + 3^{NN-1}.

Partial_Integer=zeros(MM,NN);
// Will contain, for each integer z in Integers,
// a lower aproximation of z in the basis 
// 1, 3^{1},..., 3^{NN-1}
// Partial_Integer(z,k)=s_1 + s_2*3^{1} +...+ s_k*3^{k-1}.

Future_Integer=zeros(MM,NN);
// Will contain, for each integer z in Integers,
// the image by the dynamics under the control consisting in
// extracting block in the corresponding column.

State(:,1)=pmodulo(Integers, 3^{1}) ;
// State(z,1)=s_1
Partial_Integer(:,1)=State(:,1)  ;
// Partial_Integer(z,1)=s_1
Future_Integer(:,1)= maxi(1,Integers+1-3^{1}) ;
// Dynamics (with a "maxi" because some integers in 
// Integers+1-3^{1} do not correspond to "mine profiles").
Increments(:,1)=State(:,1)-1;

for k=2:NN
  remainder= ( Integers-Partial_Integer(:,k-1) ) / 3^{k-1} ;
  // s_{k} + s_{k+1}*3 + ...
  State(:,k)=pmodulo(remainder, 3) ;
  // State(:,k)=s_{k}
  Increments(:,k)=State(:,k)-1;
  Partial_Integer(:,k)=Partial_Integer(:,k-1)+3^{k-1}*State(:,k);
// Partial_Integer(z,k)=s_1+s_2*3^{1}+...+s_k*3^{k-1}
Future_Integer(:,k)= maxi(1,Integers+3^{k-1}-3^{k}) ;
// Dynamics (with a "maxi" because some integers 
// in Integers+3^{k-1}-3^{k}
// do not correspond to "mine profiles")
end

Future_Integer(:,NN)= mini(MM,Integers+3^{NN-1}) ;
// Correction for the dynamics 
// when the last column NN is selected.
// Dynamics (with a "mini" because some integers in 
// Integers+3^{NN-1} do not correspond to "mine profiles").



// -----------------------------------------
// FROM PROFILES TO INTEGERS
// -----------------------------------------

function z=profile2integer(p)
  // p : profile as a row vector
yy= p-[ 0 p(1:($-1)) ] ;
ss=yy+1 ;
z=sum( ss .* [1 3.^{[1:(NN-1)]}] ) ;
endfunction 


// -----------------------------------------
// ADMISSIBLE INTEGERS
// -----------------------------------------

// Mine profiles are those for which 
// HH \geq p_1 \geq 0,..., HH \geq p_NN \geq 0
// that is, HH \geq y_1 + ... + y_k \geq 0 for all k
// Since, starting from the profile p(0)=(0,0,...,0)),  the 
// following extraction rule will always give "mine profiles", 
// we shall not exclude other unrealistic profiles.

Admissible=zeros(MM,NN);

Profiles= cumsum (Increments , "c" ) ; 
// Each line contains a profile, realistic or not.

adm_bottom=bool2s(Profiles<HH) ;
// A block at the bottom cannot be extracted:
// an element in adm_bottom is one if and only if
// the top block of the column is not at the bottom.

// Given a mine profile, extracting one block at the surface
// is admissible if the slope is not too high.
//
// Extracting block in column 1 is admissible 
// if and only if p_1=0.
//
// Extracting block in column j 1<j<NN) is not admissible 
// if and only if y_j=1 or y_{j+1}=-1 that is, 
// (s_j-1)=1 or (s_{j+1}-1)=-1.
// Extracting block in column j is admissible 
// if and only if s_j < 2 and s_{j+1} > 0.
//
// Extracting block in column NN is admissible 
// if and only if p_{NN}=0.

Admissible(:,1)=bool2s(Profiles(:,1)==0);
Admissible(:,NN)=bool2s(Profiles(:,NN)==0);
// Corresponds to side columns 1 and NN, for which only the 
// original top block may be extracted:
// an element in columns 1 and NN of Admissible is one 
// if and only if the pair (state,control) is admissible.

Admissible(:,2:($-1))=...
	bool2s( State(:,2:($-1))<2 & State(:,3:$)>0 ) ;
// An element in column 1<j<NN of AA is one if and only 
// s_j < 2 and s_{j+1} > 0.

Admissible=Admissible .* adm_bottom ;
// An element in column j of admissible is zero 
// if and only if 
// extracting block in column j is not admissible, 
// else it is one.


Stop_Integers=Integers( prod(1-Admissible,"c") ==1 ) ;
// Labels of states for which no decision is admissible,
// hence the decision is the retirement option


// -----------------------------------------
// INSTANTANEOUS GAIN
// -----------------------------------------

Forced_Profiles= mini(HH,  maxi(1, Profiles) ) ;
// Each line contains a profile, forced to be realistic.
// This trick is harmless and useful 
// to fill in the instantaneous gain matrix.

instant_gain=zeros(MM,NN);

for uu=1:NN
  instant_gain(:,uu)= Admissible(:,uu) .* ...
	  richness(Forced_Profiles(Future_Integer(:,uu),uu) , uu )... 
	  + (1-Admissible(:,uu)) *bottom; 
end
// When the control uu is admissible, 
// instant_gain is the richness of the top block of column uu.
// When the control uu is not admissible, instant_gain 
// has value "bottom", approximation of -infinity.

//\end{verbatim}