% 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 treating the 
%            observed proportions (relative frequencies) and total individual 
%            count (N) as constants, and randomly varying the allocation of 
%            individuals to species; see function prbcount().
%
%     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.
%             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.
%

% RE Strauss, 2/20/00

function [D,Ddiff,Dpr,E,Ediff,Epr] = diverdiff(freqs,kind,iter)
  if (nargin < 2) kind = []; end;
  if (nargin < 3) kind = []; 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
  [D,E] = diversity(freqs,kind);          % Observed index values

  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 =   zeros(C,C);
  if (get_evenness)
    Ediff = zeros(C,C);
    Epr =   zeros(C,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)                               % Randomize
    incr = 1/iter;                          % Probability increment
    f = freqs;

    for it = 1:iter                         % Iterate
      for ic = 1:C                            % For each assemblage,
        ff = freqs(:,ic);                       % isolate frequencies
        i = find(ff==0);
        if (~isempty(i))
          ff(i) = [];
        end;
        ff = ff./N(ic);                         % convert to proportions
        fp = prbcount(ff,N(ic),N(ic),1,0);      % randomize counts, given proportions
        f(:,ic) = [fp;zeros(length(i),1)]; 
      end;
      [d,e] = diversity(f,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
    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;


