mfbox_rel_hierarchical_run.m

toolbox for spm 5 for data, model free analysis

function [mfbss,params]=mfbox_rel_hierarchical_run(mfbss,params,runflag)
% run bss algorithm multiple times with different number of components
% to extract and generate tree structure
%
% Usage:
%  [mfbss,params]=mfbox_rel_hierarchical_run(mfbss,params,runflag)
%
%  mfbss   - mfbss structure (see mfbox_init)
%            name - bss algorithms
%  params  - struct with
%            threshold_type - threshold components by
%            threshold      - threshold value
%            fraction       - fraction of data to use in
%                             subsequent runs
%            compnum_type   - component number reduce function
%            factor         - factor for component_type
%            sortby         - temporalspatial
%                             temporal
%                             spatial
%                             variance
%  runflag - -1 get default parameter
%             0 interactive ask parameters
%             1 interactive ask parameters and run
%             2 run
%
% Copyright by Peter Gruber, Fabian J. Theis and Ingo R. Keck
% 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','hierarchical', ...
    struct('fraction',0.5,'factor',2,'threshold',0.01, ...
    'threshold_type','none','sortby','temporalspatial', ...
    'compnum_type','exp'));

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

if (runflag>0 && isstruct(params))
    if (runflag<3)
        prgs = mfbox_progress([],'title','BSS analysis','string', ...
            sprintf('Running BSS\nanalysis...'),'progress',[1,13]);
    end

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

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

    rparams = mfbss.params.(mfbss.name);
    if (isfield(rparams,'numComponents')==0)
        if (runflag<3), mfbox_progress(prgs,'close',[]); end
        questdlg('BSS algorithm not suitable for the analysis','Info','Ok','Ok');
        varargout{1} = mfbss;
        if (nargout>1), varargout{2} = params; end
        return
    end

    cmd_func = str2func(sprintf('mfbox_%s_run',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

    if (runflag<3)
        mfbox_progress(prgs,'string',sprintf('Running BSS\nanalysis...'),'progress',[2,13]);
    end

    sa = size(A);
    compnum = zeros(1,sa);
    n = 1;
    i = sa(2); % number of components
    while (i>1)
        switch params.compnum_type
            case 'exp'
                compnum(n) = i;
                i = max(1,floor(i/max(1.1,params.factor)));   % calculate how many times we have to do i/2 to reach 1
            case 'linear'
                compnum(n) = i;
                i = max(floor(i-params.factor));   % calculate how many times we have to do i/2 to reach 1
        end
        n = n+1;
    end
    compnum(n) = 1;
    compnum = compnum(n:-1:1);

    vAp = mfbox_calcthreshold(A,S,params.threshold_type,params.threshold, ...
        mfbss.design,reference);
    As = cell(1,n); Ws = cell(1,n); Ss = cell(1,n);
    As{n} = A(:,vAp); Ws{n} = W(vAp,:);
    Ss{n} = single(S(:,vAp));
    n = n-1;
    clear A W S

    pms = sum(mask(:));
    while (n>0)
        rparams.numComponents = compnum(n);
        if (params.fraction<1)
            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,A,W,S] = cmd_func(mfbss.X(r,:),M,rparams,max(2,runflag));
            else
                tX = mfbox_getmaskdata(mfbss.X,tg,mask);
                [f,A,W,S] = cmd_func(tX(r,:),M,rparams,max(2,runflag));
            end
            vAp = mfbox_calcthreshold(A,S,params.threshold_type,params.threshold, ...
                mfbss.design,ref);
            As{n} = A(:,vAp); Ws{n} = W(vAp,:);
            if (isempty(smask))
                Ss{n} = single(mfbss.X*Ws{n}');
            else
                Ss{n} = single(tX*Ws{n}');
            end
        else
            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);
            As{n} = A(:,vAp); Ws{n} = W(vAp,:); Ss{n} = single(S(:,vAp));
        end
        n = n-1;
        clear A W S
        if (runflag<3)
            mfbox_progress(prgs,'string',sprintf('Running BSS\nanalysis...'),'progress', ...
                [12-ceil(10*(n/length(Ss))),13]);
        end
    end

    if (runflag<3)
        mfbox_progress(prgs,'string',sprintf('Running BSS\nanalysis...'),'progress',[12,13]);
    end

    H = cell(1,length(As)-1);
    for i=1:length(As)-1
        [e,hc] = mfbox_cumindependence(As{i},As{i+1});
        [e,ho] = mfbox_overlap(Ss{i},Ss{i+1});
        v = complex(hc,ho);
        %y = mfbox_getsort(v,params.sortby);
        y = mfbox_getsort(complex(var(As{i+1}),var(Ss{i+1})),v,params.sortby);
        As{i+1} = As{i+1}(:,y);
        Ws{i+1} = Ws{i+1}(y,:);
        Ss{i+1} = Ss{i+1}(:,y);
        H{i} = v(:,y);
    end
    mfbss.extraplot = {@mfbox_plothierwin,H};

    part = cell(1,length(As)); n = 1;
    for i=1:length(As), part{i} = n:(n+size(As{i},2)-1); n = n+size(As{i},2); end
    
    Ss = cat(2,Ss{:});
    m = repmat(1:length([part{:}]),size(mfbss.X,2),1);
    v = 1./std(Ss);
    mfbss.A = cat(2,As{:})./(v(m));
    mfbss.W = cat(1,Ws{:}).*(v(m))';
    mfbss.S = struct('mask',xmask,'map',m,'dat',Ss*diag(v), ...
        'part',{part});

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