% DIVERDIFF:  Randomization test for the significance of observed differences in 
%             diversity and evenness among a set of assemblages having a common 
%             species pool.  Frequencies must be absolute (counts) rather than 
%             relative.  The possible diversity and evenness indices available 
%             are described in function diversity().
%               Significance levels (p-values) are estimated by comparing the
%             observed diversity values to a null distribution generated by
%             pooling the assemblages and randomly partitioning them into 
%             random assemblages, maintaining the original numbers of individuals
%             per assemblage.
%
%     Usage: [D,Ddiff,Dpr,E,Ediff,Epr] = diverdiff(freqs,{kind},{iter})
%
%             freqs = [S x C] matrix of absolute frequencies of S species among 
%                       C assemblages.
%             kind =  integer 1-5 (see function diversity()) indicating the 
%                       particular diversity measure to be returned [default=1].
%             iter =  optional number of randomization iterations for 
%                       significance levels [default = 0].
%             -----------------------------------------------------------------
%             D =     [1 x C] vector of diversity values.
%             Ddiff = [C x C] symmetric matrix of pairwise absolute differences 
%                       in diversity among the C assemblages.
%             Dpr =   corresponding matrix of significance levels, with ones on
%                       the diagonal.
%             E =     [1 x C] vector of evenness values.
%             Ediff = [C x C] symmetric matrix of pairwise absolute differences 
%                       in evenness among the C assemblages.
%             Epr =   corresponding matrix of significance levels, with ones on
%                       the diagonal.
%

% RE Strauss, 2/20/00
%   8/8/03 - corrected problem with default value of 'iter'.
%   8/9/03 - rewrite code for randomized probabilities.

function [D,Ddiff,Dpr,E,Ediff,Epr] = diverdiff(freqs,kind,iter)
  if (nargin < 2) kind = []; end;
  if (nargin < 3) iter = []; end;

  get_pr = 0;
  get_evenness = 0;
  if (nargout > 2) get_pr = 1; end;
  if (nargout > 3) get_evenness = 1; end;

  if (isempty(kind)) kind = 1; end;
  if (isempty(iter)) iter = 0; end;

  [S,C] = size(freqs);                    % Numbers of species and assemblages
  N = sum(freqs);                         % Total sample sizes of assemblages
  totN = sum(N);                          % Total sample size
  [D,E] = diversity(freqs,kind);          % Observed index values of assemblages

  if (iter & ~get_pr)
    disp('  DIVERDIFF warning: randomization not performed due to no output.');
    Ddiff = [];
    return;
  end;

  Ddiff = zeros(C,C);                     % Allocate output matrices
  Dpr =   eye(C);
  if (get_evenness)
    Ediff = zeros(C,C);
    Epr =   eye(C);
  else
    Ediff = [];
    Epr = [];
  end;

  for i = 1:(C-1)                         % Observed pairwise differences
    for j = 2:C                           %   (upper triangular)
      Ddiff(i,j) = abs(D(i)-D(j));
      if (get_evenness)
        Ediff(i,j) = abs(E(i)-E(j));
      end;
    end;
  end;
  
  if (iter)                               % Randomized right-tailed probabilities
    incr = 1/iter;                          % Probability increment
    rfreqs = zeros(S,C);                    % Storage for randomized frequencies
    
    grp = makegrps(1:C,N);                  % Vector of individual assemblage identities
    ind = [];
    for ic = 1:C                            % Create vector in individual species identities
      ind = [ind; makegrps(1:S,freqs(:,ic))]; % across assemblages
    end;

    for it = 1:iter                         % Iterate
      ind = ind(randperm(totN));              % Randomize species identities
      for ic = 1:C                            % For each assemblage,
        i = find(grp==ic);                      % Isolate individuals for current assemblage
        [u,f] = uniquef(ind(i),1);              % Get species frequencies
        uf = zeros(S,1);
        uf(u) = f;
        rfreqs(:,ic) = uf;                      % Plug into frequencies matrix
      end;
      [d,e] = diversity(rfreqs,kind);         % Randomized index values

      for i = 1:(C-1)                         % Randomized pairwise differences
        for j = 2:C
          ddiff = abs(d(i)-d(j));
          if (ddiff >= Ddiff(i,j))
            Dpr(i,j) = Dpr(i,j) + incr;
          end;
          if (get_evenness)
            ediff = abs(e(i)-e(j));
            if (ediff >= Ediff(i,j))
              Epr(i,j) = Epr(i,j) + incr;
            end;
          end;
        end;
      end;
    end;
  else % if (~iter)
    Dpr = [];
    Epr = [];
  end;

  if (~isempty(Ddiff))                    % Make matrices symmetric
    Ddiff = trisqmat(trilow(Ddiff'),0);
  end;
  if (~isempty(Dpr))
    Dpr = trisqmat(trilow(Dpr'),1);
  end;
  if (~isempty(Ediff))                    
    Ediff = trisqmat(trilow(Ediff'),0);
  end;
  if (~isempty(Epr))
    Epr = trisqmat(trilow(Epr'),1);
  end;

  return;


