SenT(datain, alpha)

% Sen's Trend Test with seasonallity present - Sen T
%
% Another seasonal trend test that has good power detetecting a monotonic
% trends when seasonal cycles are present. But should only be used when no
% data are missing.  When no data are missing this test is more accurate
% than the Kendall Seasonal (unless serial dependence is accounted for...
% maybe. I still need to get that part working in sktt.m function to verify).
%   - Gilbert section 17.4, page 230
%
% There is a subfunction 'rank' in this function that is used to compute 
%  ranks for all values in the dataset.  Matlab's tiedrank estimates 
%  rankings differently then required for this statistic.
%
% Syntax:
%    [T sig] = SenT(datain, alpha)
%
% inputs:
%    datain(:,1) = year  (e.g. 1999)
%    datain(:,2) = season (e.g. 1 through 12)
%    datain(:,3) = values to be analyzed
%    alpha = for two tail test (e.g. 0.05)
%
% outputs:
%    T = Sen T value
%    sig = significance using normal distribution
%
% Requirements:
%    If alpha is set to zero or not provided, significance will not be 
%    computed. To test for significance, the Statistics Toolbox is
%    required. Otherwise there are no other functions outside of Matlab
%    itself.
%
% Written by
%    Jeff Burkey
%    King County- Department of Natural Resources and Parks
%    email: [e-mail]
%    12/12/2008
function [T sig] = SenT(datain, alpha)
    [m,n] = size(datain);
    % wantplot is a flag to create a figure or not default set to no
    if exist('alpha','var') == 0
        % user didn't provide assume zero (i.e. no plot)
        alpha = 0;
        sig = nan;
    end
    
    if n ~= 3
        % there is a problem in the structure of the data
        error('ErrorTrend:SenT', 'There is a problem in the structure of the input data.\n');
    else
        % Sort by Time then by Season (or data value)
        Seasons = unique(datain(:,2));
    end
    
    NumOfSeasons = length(Seasons);
    years = unique(datain(:,1));
    nyears = length(years);
    
    % By using a base year, the user can input in data starting at any time
    % index.  It is assumed that seasons are sequential starting with the
    % number one.
    baseyear = min(datain(:,1))-1;
    % Use accumarray to put linear data into matrix nyears x nseasons
    X = accumarray([datain(:,1)-baseyear datain(:,2)],datain(:,3));
    
    % Test for missing values
    if numel(X)~= m
        % there were missing values, set outputs to nans
        error('ErrorTrend:SenT', '\nThere were missing values in the input. \nDo not use this test if there are missing values.\n');
    end
    
    % create matrix of seaonsal averages
    Xavg = repmat(mean(X),size(X,1),1);
    
    % Subtract season average from value
    X2 = X - Xavg;
    
    % linearize X2 matrix, rankings are based on entire dataset, not by
    % year and/or season.
    [r,c,v] = find(X2);
    
    % Cannot use tiedrank, see function below as substitute
    %     X3 = accumarray([r c],tiedrank(v));
    rnk = rank(v,length(v));
    
    % reconstruct ranks back into matrix
    X3 = accumarray([r c],rnk);
    
    % calculate seasonal and annual means of Ranks
    RseasonAvgs = mean(X3,1);
    RyearAvgs = mean(X3,2);
    
    % build equation for Sen's Test.
    yadd1 = (nyears + 1)/2;
    yadd2 = (nyears*NumOfSeasons + 1)/2;
    sumj = sum(sum((X3 - repmat(RseasonAvgs,nyears,1)).^2));
    sumi = sum(((years - baseyear)- yadd1).*(RyearAvgs-yadd2));
    T = sqrt(12*NumOfSeasons*NumOfSeasons/(nyears*(nyears+1)*sumj))*sumi;
    
    % Using ranks with sufficient nyears and nseasons, distribution should
    % be near normal, so test for significance using u=0 std=1
    % This ztest function requires the Statistics Toolbox
    if alpha ~= 0
        % if no alpha is provided, assume to not test.
        [h, sig] = ztest(T,0,1,alpha);
    end
    
end
% DO NOT use Matlab function tiedrank. Use this function instead.
% This function is taken out of the Gilbert rank subroutine.  It 
% could be stream lined into matrices, but for another day.
function [r nc] = rank(x, n)
    % tolerence should be user adjusted depending on the values being
    % analyzed.
    tolernce = 1.e-10;

    nc = zeros(n,1);
    r = zeros(n,1);
    
    for ii=1:n
        xi = x(ii);
        for jj=ii:n
            % if near zero assume zero, this is a floating point
            % problem commone with compilers with zero or few digits to
            % the right of the decimal.
            if abs(xi - x(jj)) < tolernce
                % assume same value
                r(jj) = r(jj) + .5;
                r(ii) = r(ii) + .5;
            else
                % not same value
                if xi > x(jj)
                    r(ii) = r(ii) + 1;
                else
                    r(jj) = r(jj) + 1;
                end
            end
        end
    end
    for ii=1:n
        ir = uint16(r(ii));
        % number of duplicate values per rank
        nc(ir) = nc(ir) + 1;
    end
end