mfbox_rel_window_run.m

toolbox for spm 5 for data, model free analysis

function [mfbss,params]=mfbox_rel_window_run(mfbss,params,runflag)
% run bss algorithm in sliding time frames
%
% Usage:
%  [mfbss,params]=mfbox_rel_window_run(mfbss,params,runflag)
%
%  mfbss   - mfbss structure (see mfbox_init)
%            name - bss algorithms
%  params  - struct with
%            threshold_type - threshold components by
%            threshold      - threshold value
%            variance       - variance threshold for concatenating windows
%            windows        - number of windows
%            overlap        - window overlap
%            spacing        - equal
%                             design
%            descomp        - design component to use for window selection
%  runflag - -1 get default parameter
%             0 interactive ask parameters
%             1 interactive ask parameters and run
%             2 run
%
% Copyright by Peter Gruber 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','window', ...
    struct('windows',10,'overlap',0.5,'spacing','equal', ...
    'threshold_type','none','threshold',0.5,'variance',0.5, ...
    'descomp',1));

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

if (runflag>0 && isstruct(params))
    if (runflag<3)
        prgs = mfbox_progress([],'title','sliding window analysis','string', ...
            sprintf('Running BSS...'),'progress',[1,2*params.windows]);
    end

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

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

    switch params.spacing
        case 'equal'
            d = [];
        case 'design'
            if (params.descomp>0)
                d = mfbss.design(:,params.descomp);
            else
                d = sum(mfbss.design,2);
            end
    end

    windows = getwindows(mfbss.timeline,params.windows,params.overlap,d);

    nw = size(windows,1);
    mw = mean(windows,2);
    [x,twindows] = min(repmat(mw.^2,1,length(windows(:)))-2*mw*windows(:)');
    twindows = reshape(twindows,size(windows));
    twindows = (twindows==repmat((1:nw)',1,size(windows,2)));
    
    cmd_func = str2func(sprintf('mfbox_%s_run',mfbss.name));
    rparams = mfbss.params.(mfbss.name);
    
    co = 0;
    cS = cell(1,nw); cA = cell(1,nw); cW = cell(1,nw);
    for w=1:nw
        if (runflag<3)
            mfbox_progress(prgs,'string',sprintf('Running time window\nanalysis...'), ...
                'progress',[w+1,2*params.windows]);
        end

        if (isempty(smask))
            [f,A,W,S] = cmd_func(mfbss.X(:,windows(w,:)),mask,rparams,max(2,runflag));
        else
            tX = mfbox_getmaskdata(mfbss.X(:,windows(w,:)),tg,mask);
            [f,A,W,S] = cmd_func(tX,mask,rparams,max(2,runflag));
        end
        if (~isempty(mfbss.design)), design = mfbss.design(:,windows(w,:));
        else design = [];
        end
        vAp = mfbox_calcthreshold(A,S,params.threshold_type,params.threshold, ...
            design,reference);
        cS{w} = single(S(:,vAp));
        cA{w} = A(:,vAp);
        cW{w} = W(vAp,:);
        co = max(co,size(cA{w},2));
        clear A W S
    end
    nbh = ((tril(ones(nw),1).*triu(ones(nw),-1)))-eye(nw);
    nbh = nbh>0;
    overlaps = cell(nw,nw,2);
    for i=1:nw
        for j=1:nw
            if (nbh(i,j))
                sect = intersect(windows(i,:),windows(j,:));
                l = any(repmat(windows(i,:), ...
                    length(sect),1)-repmat(sect',1,size(windows,2))==0);
                overlaps{i,j,1} = l;
                l = any(repmat(windows(j,:), ...
                    length(sect),1)-repmat(sect',1,size(windows,2))==0);
                overlaps{i,j,2} = l;
            end
        end
    end

    M = cell(nw,nw);
    for i=1:nw
        for j=i:nw
            if (i==j)
                % M{i,j} = diag(kurtosis(cS{i})-3);
                M{i,j} = diag(sum(cS{i}.^2>repmat(9*var(cS{i},1),size(cS{i},1),1)));
            elseif (nbh(i,j))
                % [dummy,M{i,j}] = mfbox_cumindependence(cS{i},cS{j});
                [dummy,M{i,j}] = mfbox_overlap(cS{i},cS{j},9);
                M{j,i} = M{i,j}';
            end
        end
    end

    if (runflag<3)
        mfbox_progress(prgs,'string',sprintf('Finding permutations...'), ...
            'progress',floor(params.windows*[1.5,2]));
    end

    G = mfbox_findpermutations(nbh,M,-1,2000,0.01,1);
%    G = mfbox_findpermutations(nbh,M,0.1,2000,0.01,1);
    W = cell(1,nw);
    cinds = zeros(1,nw);
    for w=1:nw
        cS{w} = cS{w}*G{w}';
        if (w>1)
            cw = (cS{w-1}'*cS{w});
            [x,y] = max(abs(cw));
            sp = sign(cw(sub2ind(size(cw),y,1:size(cw,2))));
            sp(sp==0) = 1;
            cS{w} = cS{w}*diag(sp);
            G{w} = diag(sp)*G{w};
            M{w-1,w} = G{w-1}*M{w-1,w}*G{w}';
            if (size(M{w-1,w},1)<co || size(M{w-1,w},2)<co)
                x = nan(co,co);
                x(1:size(M{w-1,w},1),1:size(M{w-1,w},2)) = M{w-1,w};
                M{w-1,w} = x;
            end
        end
        cA{w} = cA{w}*G{w}';
        cW{w} = G{w}*cW{w};
        t = size(cA{w},2);
        cinds(1:t,w) = 1:t;
    end

    for w=2:nw
        ca = cA{w-1};
        cb = cA{w};
        os = overlaps{w-1,w,1};
        ot = overlaps{w-1,w,2};
        if (sum(os)>0)
            x = mfbox_corr(ca(os,:),cb(ot,:));
        else
            x = zeros(size(ca,2),size(cb,2));
        end
        W{w-1} = complex(x,abs(M{w-1,w}));
    end
    c = 1;
    while (c<=co)
        p = find(cinds(c,:)>0,1);
        for w=p:(nw-1)
            x = cinds(c,w);
            y = cinds(c,w+1);
            l = W{w};
            if (y==0 || real(l(x,y))<params.variance)
%            if (y==0 || (real(l(x,y))<variance&&abs(imag(l(x,y)))<sqrt(variance)))
                cinds(co+1,(w+1):nw) = cinds(c,(w+1):nw);
                cinds(c,(w+1):nw) = 0;
                co = co+1;
                break;
            end
        end
        c = c+1;
    end

    for w=1:(nw-1)
        x = nan(size(W{w},1)+1,size(W{w},2)+1);
        cx = cinds(:,w); cx(cx==0) = max(cx(:))+1;
        cy = cinds(:,w+1); cy(cy==0) = max(cy(:))+1;
        x(1:size(W{w},1),1:size(W{w},2)) = W{w};
        W{w} = x(cx,cy);
    end

    rts = length(mfbss.timeline);
    mfbss.A = zeros(rts,co);
    mfbss.W = zeros(co,rts);
    seg = zeros(1,rts);
    v = zeros(rts,co);
    tc = zeros(rts,co);
    x = cumsum(max(cinds));
    x = [0,x(1:(end-1))];
    for w=1:nw
        d = windows(w,twindows(w,:));
        for i=1:size(cinds,1)
            p = cinds(i,w);
            if (p==0), continue; end
            mfbss.A(windows(w,:),i) = ...
                mfbss.A(windows(w,:),i)+cA{w}(:,p);
            mfbss.W(i,windows(w,:)) = ...
                mfbss.W(i,windows(w,:))+cW{w}(p,:);
            v(windows(w,:),i) = v(windows(w,:),i)+1;
        end
        seg(d) = w;
        tc(d,:) = x(w)+repmat(cinds(:,w)',length(d),1);
    end
    mfbss.A(v~=0) = mfbss.A(v~=0)./v(v~=0);
    mfbss.W(v'~=0) = mfbss.W(v'~=0)./v(v~=0);
    mfbss.A(isnan(mfbss.A)) = 0;
    mfbss.W(isnan(mfbss.W)) = 0;
    mfbss.S = cat(2,cS{:});
    ids = 1:size(mfbss.S,2);
    x = ~all(isnan(mfbss.S));
    ids(x) = 1:sum(x(:));
    ids(~x) = 0;
    tc(tc==0) = length(ids)+1;
    ids = [ids,sum(x(:))+1];
    mfbss.S = [mfbss.S(:,x),zeros(size(mfbss.S,1),1)];
    tc = ids(tc);

    W{nw} = {seg,windows};
    mfbss.extraplot = {@mfbox_plotwinwin,W};

    v = 1./std(mfbss.S);
    v(~isfinite(v)) = 1;
    mfbss.S = struct('mask',xmask,'map',tc,'dat', ...
        mfbss.S*diag(v));
    mfbss.A = mfbss.A./(v(tc));
    mfbss.W = mfbss.W.*(v(tc))';

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

return


function w=getwindows(l,w,o,d)

if (isempty(d))
    l = length(l);
    wl = l/(w*(1-o)+o);
    wlr = wl*(1-o);
    w = repmat(floor((0:(w-1))*wlr),ceil(wl),1)'+repmat(1:ceil(wl),w,1);
else
    d = interp1(0:(length(d)-1),d,l)';
    l = length(l);
    dn = false(size(d(:)));
    dn(2:(end-1)) = (d(2:(end-1))>d(1:(end-2)))&(d(2:(end-1))>=d(3:end));
    f = find(dn);
    if (length(f)>w)
        df = f(2:end)-f(1:(end-1));
        [x,p] = sort(df);
        dn(f(p(1:(length(f)-w))+1)) = false;
        f = find(dn);
    end
    if (min(f)>1), f = [1;f]; end
    if (max(f)<l), f = [f;l]; end
    df = f(3:end)-f(1:(end-2));
    wn = length(f)-2;
    wl = floor((1-o)*min(df)+o*l);
    w = zeros(wn,wl);
    for i=1:wn
        s = max(1,floor(f(i+1)-wl/2));
        e = s+wl-1;
        if (e>l), s = s-e+l; e = l; end
        w(i,:) = s:e;
    end
end

w = unique(w,'rows');