%% Compute Differential Scores
% Gordon Bean, January 2012
%
% Computes the differential scores.
%
% Usage
% [scores rawD] = compute_differential_scores( rawR, rawC, ... )
%
% Expects 'rawR', the reference data struct (i.e. the untreated condition),
% and 'rawC', the condition data struct (i.e. the treatment condition).
%
% Returns 'scores', a matrix containing the differential scores, and
% 'rawD', a struct containing the differential residuals, etc.
%
% Parameters
% * varMethod - string = {observed} | expErrBound | magnitudeCorrection |
%                        expErrPrior | observedBounded | none
%   The method used to compute the variance of the difference of static
%   residuals.
% 
% * smfNormFun - function handle = {@minus} | @rdivide
%   The function used to normalize out the single mutant fitness signal in
%   the plate-normalized colony sizes.
%
% * smfCalcFun - function handle = {@nanmedian} | @parzen_mode
%   The function used to estimate the single mutant fitness from the
%   plate-normalized colony sizes.
% 
% * pooledVariance - logical scalar = {1} | 0
%   Flag indicating whether the t-statistic will be calculated using a
%   pooled variance of the double mutant variance and the array mutant
%   variance.
%   

function [scores rawD params] = compute_differential_scores ...
    ( rawR, rawC, varargin )
    params = get_params( varargin{:} );
    
    if (~isfield(params, 'varmethod'))
        params.varmethod = 'observed';
    end
    
    if (~isfield(params, 'smfnormfun'))
        params.smfnormfun = @minus;
    end
    
    if (~isfield(params, 'smfcalcfun'))
        params.smfcalcfun = @nanmedian;
    end
    
    if (~isfield(params, 'pooledvariance'))
        params.pooledvariance = 1;
    end
    
    rawD.info = params;
    
    %% Compute Differentials
    rawD.size = rawC.size - rawR.size;
    rawD = add_mean_and_stdev( rawD );

    if (isfield(rawC, 'rowlabels'))
        rawD.rowlabels = rawC.rowlabels;
    end
    if (isfield(rawC, 'collabels'))
        rawD.collabels = rawC.collabels;
    end
    
    [nq na nr] = size( rawD.size );
    tmp = reshape( permute( rawD.size, [1 3 2]), [nq*nr na]);
    exp_colony_size = params.smfcalcfun(tmp,1);
    
    rawD.residuals = bsxfun( params.smfnormfun, rawD.size, ...
        exp_colony_size );
    
    obs_N = sum( ~isnan( rawD.size ), 3);
    obs_N(obs_N<2) = nan;

    cont_N = nanmedian(obs_N(:));

    %% Compute Variance
    switch lower(params.varmethod)
        case {'experrbound', 'magnitudecorrection'}
            % Expected Error Lower Bound
            if (~isfield( params, 'experr' ))
                error('You must provide ''experr''');
            end
            
            % Calculate Variance
            expvarR = nansum(params.experr(1).exp.^2,3);
            expvarC = nansum(params.experr(2).exp.^2,3);
            
            var2 = max ...
                ( (expvarR + expvarC)./(obs_N-1), nanvar(rawD.size,0,3) );
            
        case 'experrprior' % Expected Error Bayesian Prior
            if (~isfield( params, 'experr' ))
                error('You must provide ''experr''');
            end

            expvarR = nansum(params.experr(1).exp.^2,3);
            expvarC = nansum(params.experr(2).exp.^2,3);
            
            sc2 = ( expvarR + expvarC );
            dc = sum(~isnan(rawD.size),3) - 1;
            sc2 = sc2 ./ dc;
            
            var2 = (dc.*sc2 + (obs_N - 1).*nanvar(rawD.size, 0, 3)) ...
                ./ (dc + obs_N - 1);
                        
        case 'observed'
            var2 = nanvar( rawD.size, 0, 3);
            
        case 'observedbounded'
            tvar = nanvar( rawD.size, 0, 3);
            var2 = max( tvar , prctile(tvar(:), 0.01) );
            var2(var2==0) = nan;
            
        case 'none'
            var2 = ones(size(rawD.size,1), size(rawD.size,2));
            var2(var2==0) = nan;
            
        otherwise
            error('This method is not supported: %s', params.method);
        
    end
    
    if (strcmpi(params.varmethod, 'magnitudecorrection'))
        % Increase the variance for large static interactions
        
        % Compute Static Residuals
        tmp = reshape(permute(rawR.size,[1 3 2]),[nq*nr na]);
        mags.R = bsxfun(@minus, rawR.size, nanmedian(tmp,1));
        
        tmp = reshape(permute(rawC.size,[1 3 2]),[nq*nr na]);
        mags.C = bsxfun(@minus, rawC.size, nanmedian(tmp,1));
        
        mags.add = nanmean(abs(mags.R) + abs(mags.C),3);
        
        % Fit Variance
        centers = linspace( 0, prctile( mags.add(:), 100), 30);
        w = centers(2)-centers(1);
        vars = nan(length(centers),1);
        for ii = 1 : length(centers)
            list = abs( mags.add - centers(ii) ) < w;
            vars(ii) = nanmedian( sqrt(var2(list)) );
        end
        [~,eii] = min( abs(centers - prctile(mags.add(:),99.95)) );
        vars(eii:end) = vars(eii);
        
        % Calculate adjusted variances
        var2 = (sqrt(var2) ...
            ./ interp1( centers, smooth(vars), mags.add )).^2;
        
    end
    
    %% Compute Scores
    if ( params.pooledvariance )
        median_rel_error = ...
        nanmedian( in(nanstd(rawD.size,0,3) ./ rawD.meansize));
        cont_sd = max ...
            ( nanmedian(nanstd(rawD.size,0,3),1), ...
            median_rel_error * exp_colony_size );

        var1 = cont_sd.^2;

        svar = bsxfun ...
            (@plus, var2 .* (obs_N - 1), var1 * (cont_N - 1)) ...
            ./ (cont_N + obs_N - 2) .* (1/cont_N + 1./obs_N);
    else
        svar = var2;
    end
    
    scores = nanmean(rawD.residuals, 3) ./ sqrt( svar );
    
    %% FUNCTIONS -----------------------------------------------
    function raw = add_mean_and_stdev( raw )
    
        raw.meansize = nanmean( raw.size, 3 );
        raw.sdsize = nanstd( raw.size, 0, 3 );

    end

    function params = get_params( varargin )
        params = struct;

        n = length(varargin);
        if (mod(n,2) ~= 0)
            error('Uneven number of parameters and values in list.');
        end

        tmp2 = reshape(varargin, [2 n/2]);
        for jj = 1 : size(tmp2,2)
            params.(lower(tmp2{1,jj})) = tmp2{2,jj};
        end
    end
end