function pexuptail=Wilcoxrsexact(n,m,W,ranks);
% Exact upper tail p values for Wilcoxon Rank Sum statistic
% function pexuptail=Wilcoxrsexact(n,m,W,ranks);
% Borrows shamelessly from Strauss's combvals.m
% Note that Matlab's nchoosek will also generate the list
% of combinations. This program doesn't generate the full
% matrix of combinations, but calculates the test stat only.
% Input: n size of smaller group
%        m size of larger group
%        W Wilcoxon signed rank statistic
%        ranks, actual ranks of n+m items if there are ties present.
% Written by E. Gallagher, Eugene.Gallagher@umb.edu
% Help file for Strauss' combvals:
% COMBVALS: Generates the combinations of n integers taken r at a time.  The 
%           number of such combinations is given by function nc=combin().  
%     Usage: c = combvals(n,r)
%           n = number of integers (1:n) to be combined.
%           r = number to be taken at a time (0 < r <= n).
%           -------------------------------------------------------
%           c = [nc x r] matrix of combinations.

% Based on ACM Algorithm 94, J. Kurtzberg, Comm. ACM, June 1962.
% RE Strauss, 12/18/98

% An exact conditional distribution with ties follows Hollander & Wolfe p. 115
if nargin<4
    ranks=1:n+m;
    notiedr=logical(1);
else
    if length(ranks)<n+m
        error(...
            sprintf(...
            'Number of ranks (%d) doesn''t match n+m (%d)\n',...
            length(ranks),n+m));
    end
    ranks=sort(ranks);
    notiedr=logical(0);  % could do a check to see if there really are ties with ties.m
end
ranks=ranks(:);
fudranks=flipud(ranks);    
N=n+m;
r = n;
ncomb = nchoosek(N,r);  % Matlab's built-in combination function.
if W>=n*(n+m+1)-W;
    uppertail=logical(1);
else
    W=n*(n+m+1)-W;
    uppertail=logical(0);
end
if W>sum(fudranks(1:n))
    if uppertail
        error('W impossibly large')
    else
        error('W impossibly small')
    end
elseif W==sum(fudranks(1:n)) & notiedr
    if uppertail
        pexuptail=1/ncomb;
    else
        pexuptail=(ncomb-1)/ncomb;
    end
    return
end
  % Strauss's combval lists combinations in c in lexicographic
  % order, thus the critical values for sum(C) are larger than
  % observed W. We can speed up the process by using
  % Wstar=min(W,n*(m+n+1)-W) and exiting loop when Wstar fails
  % to be less than critical value
  if ncomb>1e6
      t=sprintf(...
          '%d combinations, T=%d min (1e6 combs take 1 min on p4)\n',...
          ncomb,round(ncomb/1e6)); 
      toomany=menu(t,'Stop','Continue');
      if toomany==1
          return
      end
  end
  % c = zeros(ncomb,r);     % Don't need to store values.
  Tally=0;
  j = zeros(1,r);

  for i = 1:ncomb
    b = 1;
    endflag = 0;
    while(~endflag)
      if (j(b)>=b)
        a = j(b)-b-1;
        for l = 1:b
          j(l) = l+a;
        end;
        endflag = 1;  
      else
        if (b==r)
          for b = 1:r
            j(b) = N-r-1+b;
          end;
          endflag = 1;
        end;
        b = b+1;
      end;
    end;
    % c(i,:) = N-j(r:-1:1);
    c=N-j(r:-1:1);
    if sum(ranks(c))>=W
        Tally=Tally+1;
    end
  end;
  pexuptail=Tally/ncomb;
  if ~uppertail
      pexuptail=1-pexuptail;
  end