fmri_boot_multiblock.m

绝对经典,老外制作的功能强大的matlab实现PLS_TOOBOX

function [brainlv,s,designlv,behavlv,brainscores,designscores,behavscores, ...
	lvcorrs,boot_result,behavdata,evt_list,behavdata_lst,datamatcorrs_lst,b_scores, ...
	behav_row_idx]= ...
		fmri_boot_multiblock(datamat,behavdata,evt_list, ...
		behavdata_lst, newdata_lst, num_subj_lst, ...
		num_boot,num_cond,num_subj,Clim, ...
                min_subj_per_group,is_boot_samples,boot_samples, ...
		new_num_boot,bscan)

    % Init
    %
    brainlv = [];
    s = [];
    designlv = [];
    behavlv = [];
    brainscores = [];
    designscores = [];
    behavscores = [];
    lvcorrs = [];
    boot_result = [];
    b_scores = [];
    datamatcorrs_lst = {};

    stacked_datamatcorrs = [];
    stacked_datamat = datamat;
    stacked_behavdata = behavdata;



    num_groups = length(newdata_lst);

    % keeps track of number of times a new bootstrap had to be generated
    %
    countnewtotal=0;
    num_LowVariability_behav_boots = [];
    badbeh = [];


%     boot_progress = rri_progress_ui('initialize');
%    [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond, num_boot, 0, boot_progress, ...
    [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond, num_boot, 0, ...
	min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

    if isempty(reorder)
       return;
    end;


    progress_hdl = rri_progress_ui('initialize');

    msg = 'Working on PLS ...';
    rri_progress_ui(progress_hdl, '', msg);


       k = num_cond;
       kk = length(bscan);

    % loop accross the groups, and
    % calculate datamatcorrs for each group
    %
    for i = 1:num_groups

       n = num_subj_lst(i);

       datamat = newdata_lst{i};

       rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+1/(10*num_groups));

       % compute task mean
       %
       Tdatamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);

       % compute correlation
       %
       Bdatamatcorrs = rri_corr_maps_notall(behavdata_lst{i}, datamat, n, bscan);

       rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+3/(10*num_groups));

       %  stack task and behavior - keep un-normalize data that will be
       %  used to recover the normalized one
       %
       TBdatamatcorrs = [Tdatamatcorrs; Bdatamatcorrs];

       %  stack task and behavior - normalize to unit length to reduce
       %  scaling difference
       %
       datamatcorrs = [normalize(Tdatamatcorrs,2); normalize(Bdatamatcorrs,2)];

       stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];
       datamatcorrs_lst = [datamatcorrs_lst, {Bdatamatcorrs}];

       rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+5/(10*num_groups));

    end		% for

    % Singular Value Decomposition
    %
    [r c] = size(stacked_datamatcorrs);
    if r <= c
       [brainlv, s, v] = svd(stacked_datamatcorrs',0);
    else
       [v, s, brainlv] = svd(stacked_datamatcorrs,0);
    end

    original_sal = brainlv * s;
    orig_behavlv = v * s;
    original_v = orig_behavlv;

    s = diag(s);


    rri_progress_ui(progress_hdl,'',9/10);


      %  Separate v into 2 parts: designlv and behavlv
      %
      designlv = [];
      behavlv = [];

      for g = 1:num_groups
         t = size(stacked_behavdata, 2);

         designlv = [designlv; v((g-1)*k+(g-1)*kk*t+1 : (g-1)*k+(g-1)*kk*t+k,:)];
         behavlv = [behavlv; v((g-1)*k+(g-1)*kk*t+k+1 : (g-1)*k+(g-1)*kk*t+k+kk*t,:)];
      end

      num_col = size(designlv, 2);

      % expand the num_subj for each row (cond)
      % did the samething as testvec
      %
      designscores = [];
      row_idx = [];
      last = 0;

      for g = 1:num_groups
         n = num_subj_lst(g);

         tmp = reshape(designlv((g-1)*k+1:(g-1)*k+k,:),[1, num_col*k]);
         tmp = repmat(tmp, [n, 1]);		% expand to num_subj
         tmp = reshape(tmp, [n*k, num_col]);

         designscores = [designscores; tmp];		% stack by groups

         %  take this advantage (having g & n) to get row_idx
         %
         tmp = 1:n*k;
         tmp = reshape(tmp, [n k]);
         tmp = tmp(:, bscan);

         behavdata_lst{g} = behavdata_lst{g}(tmp(:),:);

         row_idx = [row_idx ; tmp(:) + last];
         last = last + n*k;
      end

    behav_row_idx = row_idx;
    behavdata = behavdata(row_idx,:);

    % calculate behav scores
    %
    b_scores = stacked_datamat * brainlv;
    [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat(row_idx,:), ...
		stacked_behavdata(row_idx,:), ...
		brainlv, behavlv, kk, num_subj_lst);

%      lvcorrs = original_v;


    orig_corr = lvcorrs;
    [r1 c1] = size(orig_corr);
    distrib = zeros(r1, c1, num_boot+1);
    distrib(1:r1, 1:c1, 1) = orig_corr;

    rri_progress_ui(progress_hdl,'',1);

    %  Begin bootstrap loop
    %
    % reorder = rri_mkboot_order(num_cond, num_subj_lst, num_boot);

    %  move the following code up, so we don't need to run obserbation
    %  test if bootstrap can not run
    %
    % [reorder, new_num_boot] = ...
	% rri_boot_order(num_subj_lst, num_cond, num_boot);

    if new_num_boot ~= num_boot
        num_boot = new_num_boot;
        h0 = findobj(0,'tag','PermutationOptionsFigure');
        h = findobj(h0,'tag','NumBootstrapEdit');
        set(h,'string',num2str(num_boot));

        distrib = zeros(r1, c1, num_boot+1);
        distrib(1:r1, 1:c1, 1) = orig_corr;
    end

    max_subj_per_group = 8;
    if (sum(num_subj_lst <= max_subj_per_group) == num_groups)
        is_boot_samples = ones(1,num_groups);
    else
        is_boot_samples = zeros(1,num_groups);
    end

    if isempty(reorder)
       return;
    end;

    sal_sq = original_sal.^2;
    sal_sum = original_sal;
    dsal_sq = orig_behavlv.^2;
    dsal_sum = orig_behavlv;


    %  Check min% unique values for all behavior variables
    %
    num_LowVariability_behav_boots = zeros(1, size(stacked_behavdata, 2));

    for bw = 1:size(stacked_behavdata, 2)
       for p = 1:num_boot
          vv = stacked_behavdata(reorder(:,p),bw);

%          if unique(vv) <= size(stacked_behavdata, 1)*0.5
          if rri_islowvariability(vv, stacked_behavdata(:,bw))
             num_LowVariability_behav_boots(bw) = num_LowVariability_behav_boots(bw) + 1;
          end
       end
    end

    if any(num_LowVariability_behav_boots)
       disp(' ');
       disp(' ');
       disp('For at least one behavior measure, the minimum unique values of resampled behavior data does not exceed 50% of its total.');
       disp(' ');
       disp(' ');
    end


    for p = 1:num_boot

        msg = ['Working on Bootstrap:  ',num2str(p),' out of ',num2str(num_boot)];
        rri_progress_ui(progress_hdl, '', msg);
        rri_progress_ui(progress_hdl,'',p/num_boot);

        data_p = stacked_datamat(reorder(:,p),:);
        behav_p = stacked_behavdata(reorder(:,p),:);
        stacked_data = [];