function srd=srdMatrix(A,T,exact)
% Compute SRD values for matrix A
%==========================================================================
%
% Compute the scaling factor
%
nrow=size(A,1);
srd.nrow=nrow;
srd.max=srdMax(srd.nrow);
srd.scale=(100/srd.max);
%
% Get SRD values (raw and normalized)
%
fprintf('--> Compute SRD values...\n');
srd.raw=srdCompute(A,T);
srd.norm=srd.raw * (100/srd.max);
%
% Distributional information
%
fprintf('--> Compute SRD Probabilities...\n');
[srd.prob,srd.x,srd.relFreq]=srdDist(srd.nrow,srd.max,srd.raw,exact);
%    
% Rank order the SRD values
%
[~,srd.rank]=sort(srd.norm);  
%==========================================================================
end

%
%--------------------------------------------------------------------------
%
% H E L P E R   S U B - F U N C T I O N S
%
%--------------------------------------------------------------------------
%

function m=srdMax(nrow)
% Given  the number of rows, compute the maximum SRD
%==========================================================================
if rem(nrow,2)==1
    k=(nrow-1)/2;
    m=2*k*(k+1);
else
    k=nrow/2;
    m=2*k^2;
end
%==========================================================================
end

function S=srdCompute(A,T)
% Compute the SRDs, given merit matrix A and target T 
%==========================================================================
[nrow,ncol]=size(A);
[~,Tindex]=sort(T);
[~,Aindex]=sort(A);
Arank=zeros(nrow,ncol);
for i=1:nrow
    [row,~]=find(Aindex==Tindex(i));
    Arank(i,:)=row';
end  
rankIdeal=(1:nrow)';
rankDiff=abs(bsxfun(@minus,Arank,rankIdeal)); 
S=sum(rankDiff);
%==========================================================================
end

function [srdProb,x,relFreq]=srdDist(nrow,srdMax,srdRaw,exact)
% Compute SRD probability associated with each SRD value
%==========================================================================
%
% Get SRD distrubution
%
ncol=length(srdRaw);
x=0:srdMax;
idealOrder=1:nrow; 
if nrow > 10
    % Simulate 100000 SRD values 
    nSim=100000;
    S=zeros(1,nSim);
    for i=1:nSim
        randomOrder=randperm(nrow);
        rankDiff=randomOrder-idealOrder;
        S(i)=sum(abs(rankDiff));
    end
else
    % Compute exact srd
    rho=perms(1:nrow);
    S=sum(abs( bsxfun(@minus,rho,idealOrder) ),2); 
end
%
% Compute probability 
%
edges=-.5:1:srdMax+.5;
freq=histc(S,edges); freq(end)=[];
relFreq=freq/sum(freq);
srdProb=zeros(1,ncol);
if ~exact
    % Probability via normal approximation of SRD distribution
    srdProb=normcdf(srdRaw+0.5,mean(S),std(S));
else
    % Probability via SRD distribution
    P=cumsum(relFreq);
    for j=1:ncol    
        srdProb(j)=P( sum(x <= srdRaw(j)) );
    end
end
%==========================================================================
end

