cosmo corr

function c=cosmo_corr(x,y,corr_type)
% Computes correlation - faster than than matlab's "corr" for Pearson.
%
% c=comso_corr(x[,y[,corr_type]])
%
% Inputs:
%   x          PxM matrix.
%   y          PxN matrix (optional). If omitted then y=x.
%   corr_type  'Pearson' or 'Spearman' or 'Kendall' (optional). If omitted
%              then corrtype='Pearson' and the computation time is
%              significantly reduced for small matrices x and y (with
%              /tiny/ numerical imprecisions) by the use of a custom
%              implementation.
%              Using 'Kendall' required the matlab stats
%              toolbox and is currently not supported in Octave.
% Output:
%   c          MxN matrix with c(i,j)=corr(x(:,i),y(:,j),'type',corr_type).
%
% Notes:
%  - this function does not compute probability values.
%  - Using 'Kendall' for corr_type requires the matlab stats toolbox.
%
% Example:
%   % generate some pseudo-random data.
%   x=reshape(mod(2:7:100,41),5,[]);
%   y=reshape(mod(1:7:100,37),5,[]);
%   % compute builtin corr with cosmo_corr
%   % call the function first to avoid lookup delays; then measure time
%   c=corr(x,y);
%   cc=cosmo_corr(x,y);
%   % compute differences in output
%   delta=c-cc;
%   max_delta=max(abs(delta(:)));
%   fprintf('difference not greater than eps: %d\n',max_delta<=eps);
%   > difference not greater than eps: 1
%
% See also: corr
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    y_as_x=false;

    if nargin<2
        corr_type='Pearson';
        y_as_x=true;
    elseif ischar(y)
        corr_type=y;
        y_as_x=true;
    elseif nargin<3
        corr_type='Pearson';
    end

    if y_as_x
        y=x;
    end

    switch corr_type
        case 'Pearson'
            c=corr_pearson(x,y,y_as_x);

        case 'Spearman'
            x_ranks=cosmo_tiedrank(x);
            y_ranks=cosmo_tiedrank(y);

            c=corr_pearson(x_ranks,y_ranks,y_as_x);

        otherwise
            % fall-back: use Matlab's function
            % will puke if no Matlab stat toolbox
            c=corr(x,y,'type',corr_type);
    end


function c=corr_pearson(x,y,y_as_x)

    % speed-optimized version
    nx=size(x,1);
    ny=size(y,1);

    % subtract mean
    xd=bsxfun(@minus,x,sum(x,1)/nx);
    yd=bsxfun(@minus,y,sum(y,1)/ny);

    % normalization
    n=1/(size(x,1)-1);

    % standard deviation
    xs=(n*sum(xd .^ 2)).^-0.5;
    ys=(n*sum(yd .^ 2)).^-0.5;

    % compute correlations
    c=n * (xd' * yd) .* (xs' * ys);

    if y_as_x
        % ensure diagonal elements are 1
        c=(c+c')*.5;
        dc=diag(c);
        c=(c-diag(dc))+eye(numel(dc));
    end