mfbox_rel_pkmeans_run.m

toolbox for spm 5 for data, model free analysis

function [mfbss,params]=mfbox_rel_pkmeans_run(mfbss,params,runflag)
% run bss algorithm with pkmeans clustering to access reliability
%
% Usage:
%  [mfbss,params]=mfbox_rel_pkmeans_run(mfbss,params,runflag)
%
%  mfbss   - mfbss structure (see mfbox_init)
%            name - bss algorithms
%  params  - struct with
%            threshold_type - threshold components by
%            threshold      - threshold value
%            iterations     - number of repetitions
%            fraction       - fraction of the data to use in cluster
%                             steps
%  runflag - -1 get default parameter
%             0 interactive ask parameters
%             1 interactive ask parameters and run
%             2 run
%
% Copyright by Peter Gruber, Ingo R. Keck and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

error(nargchk(1,3,nargin));
error(nargchk(1,2,nargout));

if (nargin<2), params = []; end
if (nargin<3), runflag = 1; end

params = mfbox_checkparam(params,'rel','pkmeans', ...
    struct('iterations',10,'fraction',0.5,'threshold',0.01, ...
    'threshold_type','none'));

if (abs(runflag-0.5)<1)
    params = mfbox_rel_pkmeansg(mfbss,params,runflag);
end

if (runflag>0 && isstruct(params))
    if (runflag<3)
        prgs = mfbox_progress([],'title','pk-means analysis','string',sprintf('Running pkmeans\nreliability analysis...'), ...
            'progress',[1,3*params.iterations]);
    end

    [tg,mask,smask,xmask] = mfbox_getmask(mfbss.grid,mfbss.mask);

    if (~isempty(mfbss.reference)), reference = mfbss.reference(mask);
    else reference = [];
    end

    cmd_func = str2func(sprintf('mfbox_%s_run',mfbss.name));
    rparams = mfbss.params.(mfbss.name);
    if (isempty(smask))
        [f,A,W,S] = cmd_func(mfbss.X,mask,rparams,max(2,runflag));
    else
        tX = mfbox_getmaskdata(mfbss.X,tg,mask);
        [f,A,W,S] = cmd_func(tX,mask,rparams,max(2,runflag));
    end
    vAp = mfbox_calcthreshold(A,S,params.threshold_type,params.threshold, ...
        mfbss.design,reference);
    mfbss.A = A(:,vAp); mfbss.W = W(vAp,:); mfbss.S = single(S(:,vAp));

    mfbss.reliability = struct();

    if (params.iterations>0)
        sa = size(mfbss.A);
        if (params.fraction<0.9)
            transA = mfbss.A*(mfbss.A'*mfbss.A)^(-0.5);
            As = zeros([sa(2),sa(2),params.iterations+1],'single');
            As(:,:,1) = transA'*mfbss.A;
            transS = mfbss.S*(mfbss.S'*mfbss.S)^(-0.5);
            Ss = zeros([sa(2),sa(2),params.iterations+1],'single');
            Ss(:,:,1) = transS'*mfbss.S;
        else
            As = zeros([size(mfbss.A),params.iterations+1],'single');
            As(:,:,1) = mfbss.A;
        end
        pms = length(xmask);
        for i=1:params.iterations
            if (runflag<3)
                mfbox_progress(prgs,'string',sprintf('Running pkmeans\nreliability analysis...'), ...
                    'progress',params.iterations*[1,3]+[i,0]);
            end
            if (params.fraction==1)
                if (isempty(smask))
                    [f,tA,tW,tS] = cmd_func(mfbss.X,mask,rparams,max(2,runflag));
                else
                    tX = mfbox_getmaskdata(mfbss.X,tg,mask);
                    [f,tA,tW,tS] = cmd_func(tX,mask,rparams,max(2,runflag));
                end
                vAp = mfbox_calcthreshold(tA,tS,params.threshold_type,params.threshold, ...
                    mfbss.design,reference);
                tA = tA(:,vAp); tS = tS(:,vAp);
                As(:,:,i+1) = single(tA);
            elseif (params.fraction<0.9)
                r = (rand(1,pms)<params.fraction);
                M = mask; M(M) = r;
                if (~isempty(reference)), ref = reference(r);
                else ref = [];
                end
                if (isempty(smask))
                    [f,tA,tW,tS] = cmd_func(mfbss.X(r,:),M,rparams,max(2,runflag));
                else
                    tX = mfbox_getmaskdata(mfbss.X,tg,mask);
                    [f,tA,tW,tS] = cmd_func(tX(r,:),M,rparams,max(2,runflag));
                end
                vAp = mfbox_calcthreshold(tA,tS,params.threshold_type,params.threshold, ...
                    mfbss.design,ref);
                tA = tA(:,vAp); tS = tS(:,vAp);
                As(:,:,i+1) = single(transA'*tA);
                Ss(:,:,i+1) = single(transS(r,:)'*tS);
            else
                r = (rand(1,pms)<params.fraction);
                M = mask; M(M) = r;
                if (~isempty(reference)), ref = reference(r);
                else ref = [];
                end
                tX = [];
                if (isempty(smask))
                    [f,tA,tW,tS] = cmd_func(mfbss.X(r,:),M,rparams,max(2,runflag));
                else
                    tX = mfbox_getmaskdata(mfbss.X,tg,mask);
                    [f,tA,tW,tS] = cmd_func(tX(r,:),M,rparams,max(2,runflag));
                end
                vAp = mfbox_calcthreshold(tA,tS,params.threshold_type,params.threshold, ...
                    mfbss.design,ref);
                tA = tA(:,vAp);
                As(:,:,i+1) = single(tA);
            end
            clear tA tW tS
        end

        if (params.fraction<0.9)
            if (runflag<3)
                mfbox_progress(prgs,'string', ...
                    sprintf('Running pkmeans\nreliability analysis...'), ...
                    'progress',params.iterations*[2,3]);
            end

            [A,CE,CN,CL] = pkcluster(reshape(double(As),sa(2),[]),sa(2));
            mfbss.reliability.temporalpkerr = double(CE(CL(1,:))');
            mfbss.reliability.temporalpknum = CN(CL(1,:))';

            if (runflag<3)
                mfbox_progress(prgs,'string',sprintf('Running pkmeans\nreliability analysis...'), ...
                    'progress',params.iterations*[3,3]);
            end

            [A,CE,CN,CL] = pkcluster(reshape(double(Ss),sa(2),[]),sa(2));
            mfbss.reliability.spatialpkerr= double(CE(CL(1,:))');
            mfbss.reliability.spatialpknum = CN(CL(1,:))';
        else
            As = reshape(As,sa(1),[]);
            [e,d] = eig(cov(As'));
            d = diag(d); [d,x] = sort(d,'descend');
            transA = e(:,x(1:sa(2))); %TODO more reduction
            [A,CE,CN,CL] = pkcluster(transA'*As,sa(2));
            mfbss.A = transA*A;
            mfbss.W = pinv(mfbss.A);
            if (all(xmask(:)))
                mfbss.S = mfbss.X*mfbss.W';
            else
                mfbss.S = mfbss.X(xmask,:)*mfbss.W';
            end
            mfbss.reliability.pkerr = double(CE(:)');
            mfbss.reliability.pknum = double(CN(:)');
        end
    end

    if (runflag<3)
        mfbox_progress(prgs,'string',sprintf('Running pkmeans\nreliability analysis...'), ...
            'progress',params.iterations*[3,3]);
    end

    s = size(mfbss.A);
    if (s(1)>=8)
        M = zeros(s(2));
        for i=ceil(s(1)/20):ceil(s(1)/20):floor(s(1)/4)
            preA = mfbss.A(1:(end-i),:)-repmat(mean(mfbss.A(1:(end-i),:)),s(1)-i,1);
            posA = mfbss.A((1+i):end,:)-repmat(mean(mfbss.A((1+i):end,:)),s(1)-i,1);
            C = abs(diag(1./std(preA))*preA'*posA*diag(1./std(posA)))/(s(1)-i);
            rM = max(real(M),C);
            iM = imag(M);
            iM(C>=real(M)) = i;
            M = complex(rM,iM);
        end
        mfbss.extraplot = {@mfbox_plotnetworkwin,M};
    end

    m = repmat(1:size(mfbss.A,2),size(mfbss.A,1),1);
    v = 1./std(mfbss.S);
    mfbss.A = mfbss.A./(v(m));
    mfbss.W = mfbss.W.*(v(m))';
    mfbss.S = struct('mask',xmask,'map',m,'dat',mfbss.S*diag(v));

    varargout{1} = mfbss;

    if (runflag<3), mfbox_progress(prgs,'close',[]); end
end

return


function [CC,CE,CN,CL]=pkcluster(A,l)

CQ = Inf;
for i=1:10 % bad hack? in our test clustering worked the first time
    [cc,cl,cq,ce] = mfbox_pkmeans('batch',A',l);
    if (cq<CQ), CL = cl; CE = ce; CC = cc'; break; end
end
CE = accumarray(CL(:),CE(:),[l,1],@mean);
CN = accumarray(CL(:),1,[l,1]);
%CE = 1-CE;
CE = 1-(CE/max(CE(:)));
CE(CN==0) = 0;
CN = CN/max(CN(:));
%CN = CN / l;  %so CN(i)*100% gives the percentage of component found in l runs
CL = reshape(CL,[],l);