% inform_visnet.m
% Multiple neuron information analysis
% Derived from inform_visnet.c 25 Jan 2019
% This information theoretic analysis method for multiple single neuron activity is described by:
% Rolls,E.T., Treves,A. and Tovee,M.J. (1997) 
% The representational capacity of the distributed encoding of information provided by populations of neurons in the primate temporal visual cortex. 
% Experimental Brain Research 114: 149-162.
% copyright E T Rolls, Oxford Centre for Computational Neuroscience, https://www.oxcns.org
% No warranty express or implied is provided: this is research software.
% If this software is used or adapted for new published research, 
% it is requested that the original publication, which provides details, is cited.
%
% This software is made available with Rolls, E. T. (2021) Brain Computations: What and How. Oxford University Press, Oxford. https://www.oxcns.org
% Some of the coding in this program reflects the original C program
%
% Normally, this calculates the info for cells 1 to num_c, but there is a
% line below to uncomment if you wish to calculate the info for just num_c
% cells, which can be faster. Look for nc_min.
% 
% The data are read from a file 'infom_matlab.mat' the path to which may need to be set below
% This file loads num_c the number of cells, num_s the number of stimuli, and the matrix infom.
% infom has 1 row for each trial
% Each row has the stimulus number in the range 1 : num_s, the number of stimuli; followed by num_c firing rates
% The trials can be in any order. There can be different numbers of trials for each stimulus

format compact; 

warning('off', 'MATLAB:declareGlobalBeforeUse');
global decoding_method num_s pq_r  Ps nc f0 ic rc ra sd best_s; % for function Decode later in this script

max_c = 1024;   %	/* max no of cells, was 1024, then 16384 */
max_s = 100;	%	/* max no of stimuli, was 50 */

load('../VNanalysis/infom_matlab.mat'); % this loads num_c the number of cells, num_s the number of stimuli, and infom. OR
%load('infom_matlab.mat'); OR
%load('infos_matlab.mat');

%num_c = SQNET_SIZE; in VisNet
%num_s = N_GROUPS; in VisNet
fprintf('number of cells=%4d   number of stimuli=%4d\n', num_c, num_s);
if num_c > max_c
    fprintf('In inform_visnet max_c is %d and is too low. Edit inform_visnet.m.  Exiting\n', max_c);
    return
end
if num_s > max_s
    fprintf('In inform_visnet max_s is %d and is too low. Edit inform_visnet.m.  Exiting\n', max_s);
    return
end

rc = zeros(max_c, 1);		% current response */
ic = zeros(max_c, 1);		% cell code */
nc = 0;                     % cells in sample being analysed, read from file */
rcs = zeros(max_c, max_s);	% prob distr mean rates */
acs = zeros(max_c, max_s);	% prob distr nonzero rates */
ra  = zeros(max_c, max_s);	% (training set) averages */
sd  = zeros(max_c, max_s);	% (training set) variances */
ra0 = zeros(max_c,1);   	% mean avg. and var. across all stim. */
sd0 = zeros(max_c,1);
st_dev = 0;                 % global mean variance across cells */ 
f0  = zeros(max_c, max_s);  % tr. set) fract. of zero rates */
best_s = 0;                 % predicted (most probable) stimulus */
pq_r = zeros(max_s, 1);		% probabilities of each stimulus */
max_t = 50;                 % max number of trials for any one stimulus
min_t = 3;
n_tot = 0;                  % total number of trials */
x_s = zeros(max_s, 1);		% cumul. fract. of data per stimul. */
Ps = zeros(max_s, 1);		% Prob(s)    */
s0err = 0;       
read_ASCII = 1;
s_throw = 0;                % stimuli to be discarded */
niter = 0;                  % # samples of size nc */
niter_fact = 1;
pcc = 0;                    % percent correct */
decoding_method = 1;        % default (and only for Matlab) decoding method is Probability Estimation */
scramble = 0;               % scrambling of data            */
r_max = 0;  		        % max response of a cell */
s_r = 0;			        % random stimulus */
Qsq = zeros(max_s, max_s);	% frequency table */
q1 = 0;
Qq  = zeros(max_s, 1);		% frequency of each predicted s */
Qrq = zeros(max_s, max_s);	% scrambled frequency table */
q2 = 0;
Psq = zeros(max_s, max_s);	% probability table */
p1 = 0;
P2sq = zeros(max_s, max_s);	% square probabilities */
p2 = 0;
Pq  = zeros(max_s,1);		% mean assigned probability */
P2q = zeros(max_s, 1);		% mean square assigned probs */
Pr  = zeros(max_s, 1);	    % frequency of scrambled s */
Pmix = zeros(max_s, max_s); % probability  of confusion table */
I_raw = 0;                  % quantized raw info (< frequencies) */
I_boot = 0;                 % scrambled quantized info */
I_C1 = 0;		            % first correction term from data */
I_raw_sm = 0;		        % smoothed raw info (< probabilities) */
I_C1_sm = 0;		        % smoothed correction term */
I_avg_raw = zeros(max_c, 1);% average quantized info for n cells */
I_avg_boot = zeros(max_c, 1);% same but scrambled */
I_avg_C1 = zeros(max_c, 1);	% average 1st correction term */
I_avg_raw_sm = zeros(max_c, 1);% average smoothed info for n cells */
I_avg_C1_sm = zeros(max_c, 1);% average smoothed corr term */
eps = 0.000001;             % small number */
rat = zeros(max_c, 1);      % buffer for rates */
sdt = zeros(max_c, 1);      % buffer for s.devs */eps
dis = zeros(max_s, max_s);  % effective distances */
sum_p1=0.0;
sum_p2=0.0;                 % quantities to test ultrametricity SP they shold be initialized  */ 
pedantic = 0;               % SP: for a debugging  printout  */
count = 1;                  % ~ loop counter for use with stimulus array */

nc_min = 1; % this is the default, and shows the info for cells 1 to num_c, the total number of cells
nc_max = max_c;
seed = 243342;
niter_fact = 1;         % or 100; 

resfile_name = 'res';
resfile = fopen(resfile_name, 'w');
pccfile_name = 'pcc';
pccfile = fopen(pccfile_name, 'w');
outfile_name = 'infomMatlab.out';
outfile = fopen(outfile_name, 'w');

% 
% /*                  *************************************
%                     read data and set constant parameters
%                     *************************************    */ 
%                     
nc_min = num_c; % if this line is not commented, the information is calculated just for num_c cells, instead of 1 to num_c the total number of cells

nr = zeros(num_s, 1);
r = zeros(num_c, num_s, max_t); % t is the number of trials for that cell and that stimulus. r is the rate file.
n_tot = 0; % total number of trials

for lines = 1: size(infom, 1); % the number of trials
    s = infom(lines, 1); 
    for c = 1 : num_c
        rate = infom(lines, c + 1);
        t = nr(s) + 1;
        if t < (max_t + 1)
            r(c, s, t) = round(rate); % this must be an integer !!
        end
    end
    nr(s) = nr(s) + 1;
    n_tot = n_tot + 1;
end
%r
fprintf('total number of trials read = %d\n', n_tot);
% end of reading in the data -------------------------------------

count = 1; % ETR reset count here

stimulus_array = zeros(n_tot, 1);
for c = 1 : num_c
  r_max = 0.0;
  for s = 1 : num_s
     nt = nr(s);
     for t = 1: nt
        if r(c, s, t) > r_max
           r_max = r(c, s, t);
        end
     end
  end
  fprintf('max %5.1f spikes for cell %d\n',r_max,c);
  ra0(c) = 0.0;
  sd0(c) = 0.0;
  for s = 1 : num_s
      nt = nr(s);
      ra(c, s) = 0.;
      sd(c, s) = 0.;
      f0(c, s) = 0.;
      for t = 1 : nt
          ra(c, s) = ra(c, s) + r(c, s, t) / nt;
          sd(c, s) = sd(c, s) + r(c, s, t) * r(c, s, t) / (nt - 1.);
          if r(c, s, t) < eps
              f0(c, s) = f0(c, s) + 1. / nt;
          end
      end
      sd(c, s) = sd(c, s) - ra(c, s) * ra(c, s) * nt / (nt - 1.);
      ra0(c) = ra0(c) + ra(c, s) * nt / n_tot;
      sd0(c) = sd0(c) + sd(c, s) * nt / n_tot;	
      acs(c, s) = 1.;
      if ra(c, s) <= 0.2 / nt
          rcs(c, s) = 0.2 / nt;
      elseif sd(c, s) <= ra(c, s)
          rcs(c, s) = ra(c, s);
      else
          rcs(c, s) = sd(c, s) / ra(c, s) - 1. + ra(c, s);
          acs(c, s) = ra(c, s) / rcs(c, s);
      end
   end
end

x_prev = 0.;
pcc0 = 0.;
for s = 1 : num_s
    rt = nr(s) - 1;             % number of training trials */
    for s1 = 1 : num_s
        Pmix(s, s1) = 0.0;      % this is for the metric analysis */
    end
    Ps(s) = nr(s) / n_tot;
    pcc0 = pcc0 + Ps(s) * Ps(s);
    fprintf('Ps(s)=%g \n', Ps(s));
    x_s(s) = x_prev + Ps(s);
    x_prev = x_s(s);
    fprintf(outfile, ' s=%2d  nr=%2d  x_s= %g \n', s, nr(s), x_s(s));
end
fprintf('Assigned probs for each stimulus, pcc0=%4.2f\n', pcc0 * 100.0);
fprintf(resfile,'%3d %5.3f %5.3f %4.1f\n',0,0.0,0.0,100. * pcc0);
dp = 1. / n_tot;


% 
% /*                  ***********************************
%                     loop calculating info with nc cells
%                     ***********************************   
% 
%    Cross-validation procedure used (the "etr" procedure):
%    Pick at random one of the nr[s] trials for each stimulus as a test trial;
%    the remaining (rt) trials are used for "training", that is to compute
%    the average response profile across cells for each stimulus, in ra[c][s] */

if (nc_min < 1) nc_min = 1; end
if (nc_min > num_c) nc_min = num_c; end %/* Nov 2014 so that use of a large value for nc_min can be provided to calc only for num_c to save time */
if (nc_max > num_c) nc_max = num_c; end

for nc = nc_min : nc_max
%   fprintf(outfile, '\nWorking with %2d cells.\n', nc); 
    I_avg_raw(nc) = 0.;  %/* SP all this should be [nc] not [nc-1] */
    I_avg_boot(nc) = 0.;
    I_avg_C1(nc) = 0.;
    I_avg_raw_sm(nc) = 0.;
    I_avg_C1_sm(nc) = 0.;
    pcc = 0.;

%   pick niter different nc-plets */
    niter = niter_fact * (num_c - nc + 1);
    %niter = niter_fact * (num_c - nc); % test only
    for iter = 1 : niter
        %fprintf('%d th n-plet.\n', iter);
%       fprintf(outfile, '%d th n-plet.\n', iter+1);
	    icc = iter;
	    while (icc > num_c) 
            icc = icc - num_c; 
        end
%         icc
	    ic(1) = icc;
	    %for ic1 = 1 : nc - 1
        for ic1 = 1 : nc %etr test
	        ok = 0;
	        while (ok == 0)
                c = ceil(rand * num_c);
                used = 0;
                ic2 = 1;
                while (used == 0 && ic2 < ic1)
                    if (c == ic(ic2)) 
                        used = 1; 
                    end
                    ic2 = ic2 + 1;
                end
                if (used == 0)
                    ok = 1;
                    for (s = 1 : num_s) % assume s starts at 1
                        if (nr(s) <= 1)
                            ok = 0; 
                        end
                    end
                end
            end
            ic(ic1) = c;
        end
        %ic(1 : 10) % debug
        %       reset probability tables and info values */
        I_raw = 0.;
        I_boot = 0.;
	    I_C1 = 0.;
	    I_raw_sm = 0.;
	    I_C1_sm = 0.;
	    for (s = 1 : num_s)
	        Pq(s) = 0.;
	        P2q(s) = 0.;
	        Pr(s) = 0.;
	        Qq(s) = 0.;
	        for (ss = 1 : num_s)
	            Psq(s, ss) = 0.;
	            P2sq(s, ss) = 0.;
	            Qsq(s, ss) = 0.;
	            Qrq(s, ss) = 0.;
            end
        end
	 
        % run over all responses to stimulus s */
        for (s = 1 : num_s)
	        nt = nr(s);
            rt = nt - 1;	% i.e. the number of training trials */
	        for (tt = 1 : nt)
                if (pedantic == 1)
   		            fprintf(outfile,' stim %d trial %d resp',s,tt);
                end
	            % compute averages, sd's, etc. without that trial */
	            for (ic1 = 1 : nc)
		            c = ic(ic1);
		            rc(c) = r(c, s, tt);
                    if (pedantic == 1)
		                fprintf(outfile,' %4.1f ',rc(c));
                    end
		            rat(c) = ra(c, s);
		            sdt(c) = sd(c, s);
                    sd(c, s) = sd(c, s) + rat(c) * rat(c) * nt / (nt - 1.);
                    ra(c, s) = ra(c, s) - rc(c) / nt;
                    sd(c, s) = sd(c, s) - rc(c) * rc(c) / (nt - 1.);
                    ra(c, s) = ra(c, s) * nt / rt;
                    sd(c, s) = sd(c, s) * rt / (rt - 1.);
                    sd(c, s) = sd(c, s) - ra(c, s) * ra(c, s) * rt / (rt - 1.);
                    if (rc(c) < eps) f0(c, s) = f0(c, s) - 1. / nt; end
                    f0(c, s) = f0(c, s) * nt / rt;
                    acs(c, s) = 1.;
                    if (ra(c, s) <= 0.2 / rt)
                        rcs(c, s) = 0.2 / rt;
                    elseif (sd(c, s) <= ra(c, s))
                        rcs(c, s) = ra(c, s);
                    else
		                rcs(c, s) = sd(c, s) / ra(c, s) - 1. + ra(c, s);
		                acs(c, s) = ra(c, s) / rcs(c, s);
                    end
                end
                if (pedantic == 1)
		            fprintf(outfile,'\n');
		            fprintf(outfile,'Training mean rates: \n');
		            for (ic1 = 1 : nc)
			            c = ic(ic1);
			            fprintf(outfile,'cell: %d mr: ',c);
			            for (q = 1 : num_s)
			                fprintf(outfile,' %6.3f ',ra(c, q));
                        end
			            fprintf(outfile,'\n');
                    end
                end

                     
               Decode();  % Decoding function. Returns best_s

%-------------------------------------------------- End of decoding             
               

               if (pedantic ==1)
                   fprintf(outfile,'Decoded stimulus: %d \n',best_s);
                   fprintf(outfile,'***********\n');
               end

               % ~ If we're now using all cells, store best_s in stimulus_array */

               if (nc == nc_max && iter == 1)
                   % fprintf('%d %d %d\n',iter,count,best_s); */
                   stimulus_array(count) = best_s;
                   count = count + 1;
               end

               x_ran = rand;
               s_r = 1;
               while (x_ran > x_s(s_r) && s_r < num_s)
                  s_r = s_r + 1;
               end
               for (q = 1 : num_s)
                   Psq(s, q) = Psq(s, q) + dp * pq_r(q);
                   P2sq(s, q) = P2sq(s, q) + dp * pq_r(q) * pq_r(q);
                   Pq(q) = Pq(q) + dp * pq_r(q);
                   P2q(q) = P2q(q) + dp * pq_r(q) * pq_r(q);
               end
               Qsq(s, best_s) = Qsq(s, best_s) + dp;
               Qq(best_s) = Qq(best_s) + dp;
               Qrq(s_r, best_s) = Qrq(s_r, best_s) + dp;
               Pr(s_r) = Pr(s_r) + dp;
               for (ic1 = 1 : nc)
                   c = ic(ic1);
                   rc(c) = r(c, s, tt);
                   ra(c, s) = rat(c);
                   sd(c, s) = sdt(c);
                   f0(c, s) = f0(c, s) * rt / nt;
                   if (rc(c) < eps)
                       f0(c, s) = f0(c, s) + 1. / nt;
                   end
               end
           end % of tt loop

           for (ic1 = 1 : nc)
               c = ic(ic1);
               acs(c, s) = 1.;
               if (ra(c, s) <= 0.2 / nt)
                  rcs(c, s) = 0.2 / nt;
               elseif (sd(c, s) <= ra(c, s))
                  rcs(c, s) = ra(c, s);
               else
                  rcs(c, s) = sd(c, s) / ra(c, s) - 1. + ra(c, s);
                  acs(c, s) = ra(c, s) / rcs(c, s);
               end
           end
       end % end of s loop	          	    
       %       extract info values from frequencies and probabilities */   
       for (s1 = 1 : num_s)
           rt = nr(s1);
           nb = 0;
           for (s2 = 1 : num_s)
               q1 = Qsq(s1, s2);
               q2 = Qrq(s1, s2);
               p1 = Psq(s1, s2);
               if (q1 > eps)
                  I_raw = I_raw + q1 * log(q1 / (Qq(s2) * Ps(s1)));
                  nb = nb + 1;
               end
               if (q2 > eps)
                  I_boot = I_boot + q2 * log(q2 / (Qq(s2) * Pr(s1)));
               end
               if (p1 > eps)
                  I_raw_sm = I_raw_sm +  p1 * log(p1 / (Pq(s2) * Ps(s1)));
                  I_C1_sm = I_C1_sm + P2sq(s1, s2) / p1 - p1 / Ps(s1);
               end
               Pmix(s1, s2) = Pmix(s1, s2) + p1/(Ps(s1) * niter * (nc_max - nc_min + 1));
           end
           if (Qq(s1) > eps)
               I_C1 = I_C1 - 1.;
           end
           if (Pq(s1) > eps)
               I_C1_sm = I_C1_sm - (P2q(s1) / Pq(s1) - Pq(s1));
           end
           pcc = pcc + Qsq(s1, s1) / niter;

    %      determine how many decoding bins are "relevant" */
           if (nb < num_s)
               nb_x = 0.0;
               for (q = 1 : num_s)
                   qc_x = ((Qsq(s1, q) / Ps(s1)-eps)*rt + 1.)/(rt + nb);
                   if (Qsq(s1, q) > eps)
                       nb_x = nb_x + (1. - exp(log(1.-qc_x)*rt)); 
                   end
               end
               delta_N_prev = num_s*num_s;
               delta_N = (nb - nb_x)*(nb - nb_x);
               xtr = 0;
               while (delta_N < delta_N_prev && (nb+xtr) < num_s)
                   xtr = xtr + 1;
                   nb_x = 0.0;
                   gg = log(1. + .8*nb/rt)*xtr/rt;
                   %/*gg = log(1.+ nb/rt)*xtr/rt;*/
                   xxx = 0;
                   for (q = 1 : num_s)
                       if (Qsq(s1, q) > eps)
                           qc_x = (1.-gg) * ((Qsq(s1, q)/Ps(s1)) * rt + 1.)/(rt + nb);
                           nb_x = nb_x + (1. - exp(log(1.-qc_x)*rt)); 
                       elseif (xxx < xtr)
                           qc_x = gg/xtr;            
                           nb_x = nb_x + (1. - exp(log(1.-qc_x)*rt)); 
                           xxx = xxx + 1;
                       end
                   end
                   delta_N_prev = delta_N;
                   delta_N = (nb - nb_x) * (nb - nb_x);
               end
               nb = nb + xtr - 1;
               if (delta_N < delta_N_prev) 
                   nb = nb + 1;
               end
           end
           I_C1 = I_C1 + nb;
       end % end of s1 loop

       I_raw = I_raw / log(2.);
       I_raw_sm = I_raw_sm / log(2.);
       I_boot = I_boot / log(2.);
       I_C1 = I_C1 - num_s - 1.;
       I_C1 = I_C1 * dp / (2. * log(2.));
       I_C1_sm = I_C1_sm * dp / (2. * log(2.));

       fprintf(outfile, 'Cells: ');  
       for (c = 1 : nc)
           fprintf(outfile, '%d  ', ic(c));
       end
       fprintf(outfile, '\n'); 
       fprintf(outfile, 'I_ml Quant.: Raw=%5.3f Scrbld=%5.3f C1=%5.3f I=%5.3f\n', I_raw, I_boot, I_C1, (I_raw - I_C1));
       fprintf(outfile, 'I_p Smooth: Raw=%5.3f C1_sm=%5.3f I=%5.3f\n', I_raw_sm, I_C1_sm, (I_raw_sm - I_C1_sm)); 
       I_avg_raw(nc) = I_avg_raw(nc) + I_raw / niter;
       I_avg_boot(nc) = I_avg_boot(nc) + I_boot / niter;
       I_avg_C1(nc) = I_avg_C1(nc) + I_C1 / niter;
       I_avg_raw_sm(nc) = I_avg_raw_sm(nc) + I_raw_sm / niter;
       I_avg_C1_sm(nc) = I_avg_C1_sm(nc) + I_C1_sm / niter;
    end % end of iter loop

    fprintf(outfile, 'I_ml Quant.: Raw=%5.3f Scrbld=%5.3f C1=%5.3f I=%5.3f I/cell=%5.3f\n', ...
        I_avg_raw(nc), I_avg_boot(nc), I_avg_C1(nc), ...
        I_avg_raw(nc) - I_avg_C1(nc), (I_avg_raw(nc) - I_avg_C1(nc)) / nc);
    fprintf(outfile, 'I_p Smooth: Raw=%5.3f C1=%5.3f I=%5.3f I/cell=%5.3f\n', ...
        I_avg_raw_sm(nc), I_avg_C1_sm(nc), ...
        I_avg_raw_sm(nc) - I_avg_C1_sm(nc), ...
        (I_avg_raw_sm(nc) - I_avg_C1_sm(nc)) / nc);
    fprintf(resfile, '%3d %5.3f %5.3f %4.1f\n', nc, ... 
        I_avg_raw(nc) - I_avg_C1(nc), ... 
        I_avg_raw_sm(nc) - I_avg_C1_sm(nc), pcc * 100.0); 
    fprintf('%2d cells: Perc. Corr.=%4.1f;   I_ml Quant Info=%4.2f bits;   I_p Smooth Info %4.2f bits\n', ...
        nc, pcc * 100.0, I_avg_raw(nc) - I_avg_C1(nc), I_avg_raw_sm(nc) - I_avg_C1_sm(nc));
    fprintf(pccfile, '%3d %4.1f\n', nc, pcc * 100.0);

end % of nc loop


%----------------------------------------------------------

function Decode()                % We now implement Decode as a function, though it is inline code so that it can have access to all the variables

global decoding_method num_s pq_r Ps nc f0 ic rc ra sd best_s; % for function Decode

   eps = 0.001; %/*SPdec97*/ /* empiric.: this param. should be not too low */
   leps = 2.*log(1./eps);
   best_s = 0;
   p_tot = 0.;
   p_max = 0.;

   if (decoding_method == 1)	        %/* dec 1, i.e. approximate Bayesian */
                                        %/* likelihoods with Gaussian fits   */
       for (s = 1 : num_s)              %/* for each s, estimate P(r|s)P(s)  */
	       pq_r(s) = Ps(s);             %/* start by writing P(s)            */
	       for (ic1 = 1 : nc)       	%/* and multiply by the product over */
                                        %/* cells of P(rc[c]|s)              */
	           c = ic(ic1);
	           if (rc(c) < eps)	        %/* if cell c does not fire now then */
	               if (f0(c, s) > eps)	%/* if there were trials in which it */
		               pq_r(s)  = pq_r(s) * f0(c, s);	%/* also did not fire, then multiply */
                   else                 %/* by that fraction; else, by the   */
                                        %/* tail of the probab. integral     */
                       fact = phims(0., ra(c, s), sd(c, s));
                       if (fact < eps)
                           fact = eps;
                       end
                       pq_r(s) = pq_r(s) * fact; 
                   end
               else                     %/* if instead c does fire in the current vector */
	               ee = rc(c) - ra(c, s);
	               ee = ee * ee;         %/* if it fires at exactly the mean rate, then just multiply by one (i.e., do nothing)
	               if (ee > eps)       %/* if not at exactly the mean rate, then..... */
                       if (sd(c, s) < ee/leps) %/* if there was no variance in */
		                   pq_r(s) = pq_r(s) * eps; %/* train trials, then kill the test */
                       else		   %/* while if there was, do the full thing *
		                   fact = exp(-ee / (sd(c, s) * 2.0));
		                   pq_r(s) = pq_r(s) * fact;
                       end
                   end
               end
           end
	       p_tot = p_tot + pq_r(s);
	       noise_fact = 1. + eps * (rand - 0.5);
	       if (p_max < pq_r(s) * noise_fact)
	           p_max = pq_r(s) * noise_fact;
	           best_s = s;
           end
       end % of s loop 
   else
      fprintf('error - decoding method not found\n');
   end % end of the decoding method selection
                       
   for (s = 1 : num_s)
       if (p_tot == 0.0)
	       pq_r(s) = Ps(s);	%/* if they all died, equal honor to all */
	       best_s = ceil(num_s * rand);
       else
	       pq_r(s) = pq_r(s) / p_tot;	%/* finally normalize to get P(s|r)  */
       end
   end
end %  of Decode function                         

