mfbox_tool_compare_run.m

toolbox for spm 5 for data, model free analysis

function ret=mfbox_tool_compare_run(varargin)
% compare components from different bss runs
%
% Usage:
%  mfbss_compare=mfbox_tool_compare_run(mfbss_compare,[keyword,value[...]])
%  mfbss_compare=mfbox_tool_compare_run(mfbss[,mfbss[...]],[keyword,value[...]])
%
%  mfbss_compare - mfbss compare structure
%  mfbss         - mfbss structure (see mfbox_init)
%  keyword       - valid keywords and value
%                  plot_comb
%                  plot_diff
%                  plot1
%                  plot2
%                  spat_auto_group
%                  temp_auto_group
%                  spatiotemp_plot
%                  spatial_match
%                  temporal_match
%                  add_data
%                  remove_data
%                  dataselect1
%                  dataselect2
%                  run_ica
%                  save_comparedata
%                  load_comparedata
%
% 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

handles.filename = {};
handles.mfbss = {};
handles.info = {};
handles.private = {};
handles.spatial_cross_data = []; handles.temporal_cross_data = [];
handles.grouping = {};
handles.components = struct('timeline',[]);

handles.frontend = struct();
handles.newfrontend = struct();
handles.newimages = struct();
handles.newplots = struct();
handles.buttons = struct();
handles.interactive = 1;
handles.noninteractive = {};
ret = struct();

if (nargin>0)
    num = 1;
    p = 0;
    action = '';
    for i=1:nargin
        if (p==0 && isstruct(varargin{i}))
            if (isfield(varargin{i},'mfbss'))
                comparehandles = varargin{i};
                fields = {'mfbss','filename','private','info','spatial_cross_data', ...
                    'temporal_cross_data','components','components_masks','grouping'};
                for j=1:length(fields)
                    if (isfield(comparehandles,fields{j}))
                        handles.(fields{j}) = comparehandles.(fields{j});
                    end
                end
            else
                [handles.mfbss{num},handles.filename{num},handles.private{num},handles.info{num}] = ...
                    load_data(varargin{num},'',struct(),handles.interactive);
                num = num+1;
            end
        else
            if (p==0), p = i; end
            if (mod(i-p,2)==0)
                action = varargin{i};
            else
                param = varargin{i};
                switch action
                    otherwise
                        warning('unknown parameter %s',action);
                end
            end
        end
    end
    handles = recalcinfo(handles,1,num);
end

if (handles.interactive)
    handles.runflag = 0;
    task = 'start';
    while (handles.runflag<2)
        [fig,task,buttons,data] = mfbox_tool_compareg(handles.buttons, ...
            struct('info',{handles.info},'grouping',{handles.grouping}, ...
            'mfbsscomp',{cellfun(@(x) size(x.A,2),handles.mfbss)}, ...
            'components_timeline',cat(2,handles.components.timeline), ...
            'frontend',handles.newfrontend,'plots',handles.newplots, ...
            'images', handles.newimages),task);
        handles.newfrontend = struct();
        handles.newplots = struct();
        handles.newimages = struct();
        if (handles.runflag==0)
            handles.runflag = 1;
            handles.frontend = data.frontend;
            handles.grouping = data.grouping;
            handles = reset_views(handles);
        elseif (strcmp(task,'stop'))
            handles.runflag = 2;
        else
            handles.frontend = data.frontend;
            handles.grouping = data.grouping;
            button = false;
            b = fieldnames(buttons);
            for i=1:length(b)
                if (buttons.(b{i})==1)
                    try
                        cmd_func = str2func(b{i});
                        handles = cmd_func(handles);
                    catch
                        warning('something bad happened while executing %s (%s)', ...
                            b{i},lasterror());
                    end
                    button = true;
                end
            end
            if (button)
                handles = reset_views(handles);
            else
                handles = redraw_all(handles);
            end
        end
    end
else
    for i=1:2:length(noninteractive(:))
        try
            cmd_func = str2func(noninteractive{i});
            handles = cmd_func(handles,noninteractive{i+1});
        catch
            warning('something bad happened while executing %s (%s)', ...
                noninteractive{i},lasterror());
        end
    end
end

if (nargout>0)
    fields = {'mfbss','filename','private','info','spatial_cross_data', ...
        'temporal_cross_data','components','components_masks','grouping'};
    for i=1:length(fields)
        if (isfield(handles,fields{i}))
            ret.(fields{i}) = handles.(fields{i});
        end
    end
end

return


% --- redraw all diagrams
function handles=redraw_all(handles)

curgroup = handles.frontend.groups;
if (length(handles.mfbss)>=2)
    names = mknames(length(handles.mfbss),'Data');
    sel1 = handles.frontend.dataselect1;
    sel2 = handles.frontend.dataselect2;
    num = length(handles.mfbss);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'dataselect1','String',names(sel2~=(1:num)), ...
        'Value',sel1-(sel2<sel1));
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'dataselect2','String',names(sel1~=(1:num)), ...
        'Value',sel2-(sel1<sel2));
    agroup = handles.mfbss{sel1}.S.tag;
    bgroup = handles.mfbss{sel2}.S.tag;
elseif (length(handles.mfbss)==1)
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'dataselect1','String',{'Data1'},'Value',1);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'dataselect2','String',{''},'Value',1);
    sel1 = 1;
    sel2 = 0;
    agroup = handles.mfbss{sel1}.S.tag;
    bgroup = {};
else
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'components1','String',{' '},'Value',1);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'components2','String',{' '},'Value',1);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'groups','String',{' '},'Value',1);
    return;
end
g = handles.grouping;
handles.newfrontend = addfrontend(handles.newfrontend, ...
    'components1','String',agroup);
handles.newfrontend = addfrontend(handles.newfrontend, ...
    'components2','String',bgroup);

if (size(g,2)==0), groups = {' '};
else groups = mknamespat(1:size(g,2),'Group %d');
end

curgroup = min(length(groups),max(1,curgroup));
handles.newfrontend = addfrontend(handles.newfrontend, ...
    'groups','String',groups,'Value',curgroup);
if (~isempty(g))
    if (~isempty(agroup))
        handles.newfrontend = addfrontend(handles.newfrontend, ...
            'components1','Value',g{sel1,curgroup});
    end
    if (~isempty(bgroup))
        handles.newfrontend = addfrontend(handles.newfrontend, ...
            'components2','Value',g{sel2,curgroup});
    end
end

if (size(g,2)>=curgroup)
    ma = -Inf; mi = Inf;
    plotdata = {};
    for i=1:size(g,1)
        if ((i~=sel1)&&(i~=sel2))
            plotdata{end+1} = {handles.components(i).timeline, ...
                handles.components(i).tcs, ...
                'Color',[0.6;0.6;0.6],'LineWidth',6,'LineStyle','--'};
            t = handles.components(i).tcs(:,g{i,curgroup});
            if (~isempty(t))
                plotdata{end+1} = {handles.components(i).timeline,t, ...
                    'Color',[0;0;0],'LineWidth',2};
                mi = min(mi,min(t(:))); ma = max(ma,max(t(:)));
            end
        end
    end
    t = []; v = [];
    if (sel1>0)
        plotdata{end+1} = {handles.components(sel1).timeline, ...
            handles.components(sel1).tcs, ...
            'Color',[0.8;0.8;0.6],'LineWidth',3,'LineStyle','-'};
        t = handles.components(sel1).tcs(:,g{sel1,curgroup});
    end
    if (sel2>0)
        plotdata{end+1} = {handles.components(sel2).timeline, ...
            handles.components(sel2).tcs, ...
            'Color',[0.6;0.8;0.8],'LineWidth',3,'LineStyle','-'};
        v = handles.components(sel2).tcs(:,g{sel2,curgroup});
    end
    if (~isempty(t))
        plotdata{end+1} = {handles.components(sel1).timeline,t, ...
            'Color',[1;0.2;0.2],'LineWidth',1.5};
        mi = min(mi,min(t(:))); ma = max(ma,max(t(:)));
    end
    if (~isempty(v))
        plotdata{end+1} = {handles.components(sel2).timeline,v, ...
            'Color',[0.2;0.2;1],'LineWidth',1.5};
        mi = min(mi,min(v(:))); ma = max(ma,max(v(:)));
    end
    axesdata = {'YLim',[mi,ma]+[-1,1]*(ma-mi)/2};
    handles.newplots.timecourse = struct('plots',{plotdata},'axes',{axesdata});
    xmin = min(cat(2,handles.components.timeline));
    xmax = max(cat(2,handles.components.timeline));
    sv = handles.frontend.position;
    zf = handles.frontend.zoomfact;
    inc = (xmax-xmin)/zf;
    xmax = xmax-inc;
    sv = min(max(sv,xmin),xmax);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'position','Min',xmin,'Max',xmax,'Value',sv);
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'timecourse','XLim',[sv,sv+inc]);
    l = [0 cumsum(cellfun(@(x) size(x.A,2),handles.mfbss))];
    h = g;
    for i=1:size(h,1)
        for j=1:size(h,2), h{i,j} = l(i)+h{i,j}; end
    end
    [dummy,sinv] = sort(cat(2,h{:}));
    [dummy,s] = sort(sinv);
    os = handles.frontend.orderselect;
    if (os==1)
        ss = real(handles.spatial_cross_data(s,s));
        tt = real(handles.temporal_cross_data(s,s));
        handles.newimages.allcorr = {reshape(reshape(tril(ss),[],1)*[0,0,1]+ ...
            reshape(triu(tt),[],1)*[1,0,0],l(end),l(end),3)};
    elseif (os==2)
        ss = imag(handles.spatial_cross_data(s,s));
        tt = imag(handles.temporal_cross_data(s,s));
        handles.newimages.allcorr = {reshape(reshape(tril(ss),[],1)*[0,0,1]+ ...
            reshape(triu(tt),[],1)*[1,0,0],l(end),l(end),3)};
    end
end

if ((sel1==0)||(sel2==0)), return; end

l = [0 cumsum(cellfun(@(x) size(x.A,2),handles.mfbss))];
mark = zeros(l(sel1+1)-l(sel1),l(sel2+1)-l(sel2));
mark(g{sel1,curgroup},:) = mark(g{sel1,curgroup},:)+1;
mark(:,g{sel2,curgroup}) = mark(:,g{sel2,curgroup},:)+1;
mark(mark>0) = 0.5;
if (~isempty(handles.spatial_cross_data))
    os = handles.frontend.orderselect;
    if (os==1)
        t = reshape(real(handles.spatial_cross_data((l(sel1)+1):l(sel1+1), ...
            (l(sel2)+1):l(sel2+1))),[],1);
    elseif (os==2)
        t = reshape(imag(handles.spatial_cross_data((l(sel1)+1):l(sel1+1), ...
            (l(sel2)+1):l(sel2+1))),[],1);
    else
        t = zeros(1);
    end
    t = t/max(t(:));
    t = reshape(t*[0,0,1]+mark(:)*[0,0.5,0],[size(mark),3]);
    handles.newimages.spatial_cross = {t};
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'spatial_cross','YDir','reverse');
end
if (~isempty(handles.temporal_cross_data))
    if (os==1)
        t = reshape(real(handles.temporal_cross_data((l(sel1)+1):l(sel1+1), ...
            (l(sel2)+1):l(sel2+1))),[],1);
    elseif (os==2)
        t = reshape(imag(handles.temporal_cross_data((l(sel1)+1):l(sel1+1), ...
            (l(sel2)+1):l(sel2+1))),[],1);
    else
        t = zeros(1);
    end
    t = t/max(t(:));
    t = reshape(t*[1,0,0]+mark(:)*[0,0.5,0],[size(mark),3]);
    handles.newimages.temporal_cross = {t};
    handles.newfrontend = addfrontend(handles.newfrontend, ...
        'temporal_cross','YDir','reverse');
end


% --- redraw all diagrams
function handles=reset_views(handles)

buttons = struct('plot_comb',0,'plot_diff',0,'plot1',0,'plot2',0, ...
    'spat_auto_group',0,'temp_auto_group',0,'spatiotemp_plot',0, ...
    'spatial_match',0,'temporal_match',0,'add_data',1, ...
    'remove_data',0,'dataselect1',0,'dataselect2',0,'run_ica',0, ...
    'save_comparedata',0,'load_comparedata',1);

if (length(handles.mfbss)>=2)
    varthresh_opts = {'>1 std, 0 centered','>3 std, 0 centered', ...
        '>1 std, mean centered','>3 std, mean centered'};
    if (isempty(cat(2,handles.components.timeline)))
        l = [0 cumsum(cellfun(@(x) size(x.A,2),handles.mfbss))];
        handles.spatial_cross_data = zeros(l(end));
        handles.temporal_cross_data = zeros(l(end));
        handles.components_masks = struct();
        handles.components = struct('timeline',[]);
        i = 1; j = 2; dat = {};
        while((i<=length(handles.mfbss))&&(j<=length(handles.mfbss)))
            [timeline,tcs1,tcs2] = gettcs( ...
                handles.mfbss{i}.A,handles.mfbss{i}.timeline, ...
                handles.mfbss{j}.A,handles.mfbss{j}.timeline);
            Sa = zeros(size(handles.mfbss{i}.S.dat,1),size(handles.mfbss{i}.A,2));
            for ii=1:size(Sa,2)
                Sa(:,ii) = sum(handles.mfbss{i}.S.dat(:,unique(handles.mfbss{i}.S.map(:,ii))),2);
            end
            Sb = zeros(size(handles.mfbss{j}.S.dat,1),size(handles.mfbss{j}.A,2));
            for jj=1:size(Sb,2)
                Sb(:,jj) = sum(handles.mfbss{j}.S.dat(:,unique(handles.mfbss{j}.S.map(:,jj))),2);
            end
            [mask,dat1,dat2,handles.private{i},handles.private{j}] = getdat( ...
                handles.private{i},handles.mfbss{i}.grid(:,handles.mfbss{i}.S.mask), ...
                Sa,handles.private{j}, ...
                handles.mfbss{j}.grid(:,handles.mfbss{j}.S.mask),Sb);
            [sc,tc] = update_cross(timeline,tcs1,tcs2,mask,dat1,dat2);
            handles.spatial_cross_data((l(i)+1):l(i+1),(l(j)+1):l(j+1)) = sc;
            handles.temporal_cross_data((l(i)+1):l(i+1),(l(j)+1):l(j+1)) = tc;
            handles.spatial_cross_data((l(j)+1):l(j+1),(l(i)+1):l(i+1)) = sc.';
            handles.temporal_cross_data((l(j)+1):l(j+1),(l(i)+1):l(i+1)) = tc.';
            for k=1:(i-2)
                [timeline,tcs1,tcs2] = gettcs( ...
                    handles.mfbss{k}.A,handles.mfbss{k}.timeline, ...
                    handles.mfbss{i}.A,handles.mfbss{i}.timeline);
                [sc,tc] = update_cross(timeline,tcs1,tcs2,mas{k}+2*bitand(mask,1),dat{k},dat1);
                handles.spatial_cross_data((l(k)+1):l(k+1),(l(i)+1):l(i+1)) = sc;
                handles.temporal_cross_data((l(k)+1):l(k+1),(l(i)+1):l(i+1)) = tc;
                handles.spatial_cross_data((l(i)+1):l(i+1),(l(k)+1):l(k+1)) = sc.';
                handles.temporal_cross_data((l(i)+1):l(i+1),(l(k)+1):l(k+1)) = tc.';
            end
            for k=1:(j-2)
                [timeline,tcs1,tcs2] = gettcs( ...
                    handles.mfbss{k}.A,handles.mfbss{k}.timeline, ...
                    handles.mfbss{j}.A,handles.mfbss{j}.timeline);
                [sc,tc] = update_cross(timeline,tcs1,tcs2,mas{k}+bitand(mask,2),dat{k},dat2);
                handles.spatial_cross_data((l(k)+1):l(k+1),(l(j)+1):l(j+1)) = sc;
                handles.temporal_cross_data((l(k)+1):l(k+1),(l(j)+1):l(j+1)) = tc;
                handles.spatial_cross_data((l(j)+1):l(j+1),(l(k)+1):l(k+1)) = sc.';
                handles.temporal_cross_data((l(j)+1):l(j+1),(l(k)+1):l(k+1)) = tc.';
            end
            [timeline,dummy,tcs1] = gettcs([],[], ...
                handles.mfbss{i}.A,handles.mfbss{i}.timeline);
            [timeline,dummy,tcs2] = gettcs([],[], ...
                handles.mfbss{j}.A,handles.mfbss{j}.timeline);
            [sc,tc] = update_cross(timeline,tcs1,tcs1,3*bitand(mask,1),dat1,dat1);
            handles.spatial_cross_data((l(i)+1):l(i+1),(l(i)+1):l(i+1)) = sc;
            handles.temporal_cross_data((l(i)+1):l(i+1),(l(i)+1):l(i+1)) = tc;
            handles.components(i).mask = mask==1|mask==3;
            m = mean(dat1); s = std(dat1);
            d = zeros(size(dat1),'uint8');
            for k=1:length(m)
                d(:,k) = (abs(dat1(:,k))>=s(k))+2*(abs(dat1(:,k))>=3*s(k))+ ...
                    4*(abs(dat1(:,k)-m(k))>=s(k))+8*(abs(dat1(:,k)-m(k))>=3*s(k));
            end
            handles.components(i).dat = d;
            handles.components(i).timeline = timeline;
            handles.components(i).tcs = tcs1;
            [sc,tc] = update_cross(timeline,tcs2,tcs2,3*bitand(mask,2)/2,dat2,dat2);
            handles.spatial_cross_data((l(j)+1):l(j+1),(l(j)+1):l(j+1)) = sc;
            handles.temporal_cross_data((l(j)+1):l(j+1),(l(j)+1):l(j+1)) = tc;
            handles.components(j).mask = mask==2|mask==3;
            m = mean(dat2); s = std(dat2);
            d = zeros(size(dat2),'uint8');
            for k=1:length(m)
                d(:,k) = (abs(dat2(:,k))>=s(k))+2*(abs(dat2(:,k))>=3*s(k))+ ...
                    4*(abs(dat2(:,k)-m(k))>=s(k))+8*(abs(dat2(:,k)-m(k))>=3*s(k));
            end