fmri_taskpls_norotate.m

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

function [brainlv, s, designlv, b_scores, d_scores, lvintercorrs, design, ...
           perm_result, boot_result, dev_evt_list] = ...
		fmri_taskpls_norotate(st_datamat, design, num_conditions,...
			evt_list, num_boot, num_perm, subj_group, ...
         min_subj_per_group,is_boot_samples,boot_samples,new_num_boot)

  progress_hdl = ShowProgress('initialize','Working on PLS ... ');

  msg = sprintf('Working on PLS ... ');
  ShowProgress(progress_hdl,msg)
  if ~isempty(progress_hdl)
     setappdata(progress_hdl,'ProgressScale',1/(num_perm+1)*0.8);
     setappdata(progress_hdl,'ProgressStart',0.2);
     ShowProgress(progress_hdl,0)
  end;


  if isempty(subj_group),
     GROUP_ANALYSIS = 0;			% for nongroup analysis
  else
     GROUP_ANALYSIS = 1;
  end;

  if (GROUP_ANALYSIS == 0),			
%     dev_data = st_datamat - ones(size(st_datamat,1),1)*mean(st_datamat,1);
%     dev_evt_list = evt_list;
     [dev_data dev_design] = ...
	gen_dev_data(st_datamat,num_conditions,evt_list,design);
     dev_evt_list = evt_list;
  else
     [dev_data dev_design] = ...
	gen_grp_dev_data(st_datamat,num_conditions,evt_list,subj_group,design);
%     dev_evt_list = repmat([1:num_conditions],1,length(subj_group));
     dev_evt_list = evt_list;
  end;


  %  instead of SVD
  crossblock = normalize(dev_design)'*dev_data;
  brainlv = crossblock';
  s = sqrt(sum(crossblock.^2, 2));
  designlv = dev_design;
  normalized_brainlv = normalize(brainlv);
  lvintercorrs = normalized_brainlv'*normalized_brainlv;


  if (GROUP_ANALYSIS == 0),			
     d_scores = designlv(evt_list,:);
  else
     d_scores = designlv;

     %  need to expand d_score and reorder with evt_list grp by grp
     %
     d_scores = expand_d_scores(d_scores, num_conditions, evt_list, subj_group);
  end;
  b_scores = st_datamat * brainlv;


  perm_result = [];
  boot_result = [];


  if num_perm > 0

    %  Begin permutation loop
    %
    ShowProgress(progress_hdl,1);
    msg = sprintf('Computing permutation orders ...');
    ShowProgress(progress_hdl,msg);

    % generate the permutation orders
    %
    if (GROUP_ANALYSIS == 0),			
       num_repetitions = length(evt_list) / num_conditions;
       perm_order = rri_perm_order(num_repetitions,num_conditions,num_perm);
    else
       perm_order = rri_perm_order(subj_group,num_conditions,num_perm);
    end;

    permcount = zeros(size(s));

    for k=1:num_perm,

      msg = sprintf('Working on permutation ... %d out of %d',k,num_perm);
      ShowProgress(progress_hdl,msg);

      new_order = perm_order(:,k);
      if (GROUP_ANALYSIS == 0),			% for nongroup analysis 
         [data_p design_p] = gen_dev_data(st_datamat(new_order,:),num_conditions, ...
					evt_list,design);
      else
         [data_p design_p] = gen_grp_dev_data(st_datamat(new_order,:),num_conditions, ...
					evt_list,subj_group,design);
      end;

      crossblock = normalize(design_p)'*data_p;
      s_perm = sqrt(sum(crossblock.^2, 2));
      permcount = permcount + (s_perm >= s);;

      ShowProgress(progress_hdl,k+1);

    end; 
    %
    %-- permutation loop ------------------------------------------------- 

    perm_result.s_prob = permcount / num_perm;
    perm_result.num_perm = num_perm;
    perm_result.permsamp = perm_order;
    perm_result.sp = permcount;
    perm_result.dp = [];
    perm_result.designlv_prob = [];

  end					% end permutation loop


  if num_boot > 0

    %  Begin bootstrap loop
    %
    max_subj_per_group = 8;

    if isempty(subj_group),
       num_groups = 0;
       GROUP_ANALYSIS = 0;			% for nongroup analysis

       %  moved from below
       %
       num_rep = length(evt_list) / num_conditions;
 
       % make sure the datamat are blocked by subjects
       [sort_evt,sort_idx] = sort(evt_list);
       new_row_list = reshape(sort_idx,num_rep,num_conditions)';
       new_evt_list = evt_list(new_row_list);

%       boot_progress = rri_progress_ui('initialize');
%       [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,0,boot_progress, ...
       [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,0, ...
           min_subj_per_group,is_boot_samples,boot_samples,new_num_boot); 

       if num_rep <= max_subj_per_group
          is_boot_samples = 1;
       else
          is_boot_samples = 0;
       end

    else
       num_groups = length(subj_group);
       GROUP_ANALYSIS = 1;

       %  moved from below
       %
%       boot_progress = rri_progress_ui('initialize');
%       [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,0,boot_progress, ...
       [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,0, ...
           min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

       if (sum(subj_group <= length(subj_group)) == num_groups)
          is_boot_samples = 1;
       else
          is_boot_samples = 0;
       end

    end;

    if isempty(boot_order)
       boot_result = [];
       return;
    end

    ShowProgress(progress_hdl,1);
    msg = sprintf('Compute bootstrap resampling orders ...');
    ShowProgress(progress_hdl,msg);

    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));

       if ~isempty(progress_hdl)
          setappdata(progress_hdl,'ProgressScale',1/(num_boot+1)*0.8);
       end
    end

    if isempty(boot_order)
       boot_result = [];
       return;
    end


    brainlv_sq = brainlv.^2;
    brainlv_sum = brainlv;


    for k=1:num_boot,

      msg = sprintf('Working on bootstrap ... %d out of %d',k,num_boot);
      ShowProgress(progress_hdl,msg);

      new_order = boot_order(:,k);

      if (GROUP_ANALYSIS == 0),			% for nongroup analysis
         new_row_order = new_row_list(:,new_order);
         [data_p design_p] = gen_dev_data(st_datamat(new_row_order,:), ...
			num_conditions, new_evt_list, design);
      else
         [data_p design_p]= gen_grp_dev_data(st_datamat(new_order,:), ...
			num_conditions, evt_list, subj_group, design);
      end;

      crossblock = normalize(design_p)'*data_p;
      brainlv_sq = brainlv_sq + (crossblock.^2)';
      brainlv_sum = brainlv_sum + crossblock';
    
      ShowProgress(progress_hdl,k+1);

    end; 
    %
    %-- bootstrap loop ------------------------------------------------- 

    boot_result.num_boot = num_boot;
    boot_result.bootsamp = boot_order;
    brainlv_sum2 = (brainlv_sum.^2) / (num_boot + 1);
    boot_result.brainlv_se=sqrt((brainlv_sq - brainlv_sum2)/num_boot);

    %  check for zero standard errors - replace with ones
    %
    test_zeros=find(boot_result.brainlv_se<=0);

    boot_result.zero_brain_se = test_zeros;

    if ~isempty(test_zeros);
        boot_result.brainlv_se(test_zeros)=1;
    end

    boot_result.compare = brainlv ./ boot_result.brainlv_se;

    %  for zero standard errors - replace bootstrap ratios with zero
    %  since the ratio makes no sense anyway
    %
    if ~isempty(test_zeros);
         boot_result.compare(test_zeros)=0;
    end

  end					% end bootstrap loop

  return;						% fmri_taskls_norotate


%-------------------------------------------------------------------------
function  [data, dev_design] = ...
	gen_grp_dev_data(st_datamat, num_conditions, evt_list, subj_group, design)
%
%  compute the average of st_datamat within the same group
%
  num_group = length(subj_group);
  g_end_idx = 0;
  data = [];
  dev_design = [];

  for g = 1:num_group,

     g_start_idx = g_end_idx + 1;
     g_end_idx = g_start_idx+subj_group(g)*num_conditions-1;
     g_range = [g_start_idx:g_end_idx];

     %  store the row indices for each task 
     %
     g_evt_list = zeros(1,length(evt_list));

     g_evt_list(g_range) = evt_list(g_range);

     task_idx = cell(1,num_conditions);
     for i=1:num_conditions,
        task_idx{i} = find(g_evt_list==i);
     end; 

     %  compute the mean data of each condition for the group 
     %
     if strcmpi(class(st_datamat),'single')
        mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
     else
        mean_datamat = zeros(num_conditions,size(st_datamat,2));
     end

     for i=1:num_conditions,
        mean_datamat(i,:) =  mean(st_datamat(task_idx{i},:),1);
     end;
     grp_datamat = mean_datamat; % - ones(num_conditions,1)*mean(mean_datamat,1);
     data = [data; grp_datamat];

     span = [((g-1)*num_conditions+1) : g*num_conditions];
     dev_design = [dev_design; design(span,:)];
  end;

  return; 					% gen_grp_dev_data


%-------------------------------------------------------------------------
function  [data, design] = ...
	gen_dev_data(st_datamat, num_conditions, evt_list, design)

  task_idx = cell(1,num_conditions);
  for i=1:num_conditions,
     task_idx{i} = find(evt_list==i);
  end; 

  %  compute the mean data of each condition for the group 
  %
  if strcmpi(class(st_datamat),'single')
     mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
  else
     mean_datamat = zeros(num_conditions,size(st_datamat,2));
  end

  for i=1:num_conditions,
     mean_datamat(i,:) =  mean(st_datamat(task_idx{i},:),1);
  end;

  data = mean_datamat; % - ones(num_conditions,1)*mean(st_datamat,1);

  return; 					% gen_dev_data


%-------------------------------------------------------------------------
function hdl = ShowProgress(progress_hdl,info)

  %  'initialize' - return progress handle if any
  %
  if ischar(progress_hdl) & strcmp(lower(progress_hdl),'initialize'),
     if ~isempty(gcf) & isequal(get(gcf,'Tag'),'ProgressFigure'),
         hdl = gcf;
         if ~isempty(info)
             set(hdl,'Name',info);
         end;
     else
         hdl = [];
     end;
     return;
  end;


  if ~isempty(progress_hdl)
     if ischar(info)
         rri_progress_status(progress_hdl,'Show_message',info);
     else
         rri_progress_status(progress_hdl,'Update_bar',info);
     end;
     return;
  end;

  if ischar(info),
     disp(info)
  end;

  return;					% ShowProgress


%-------------------------------------------------------------------------
function new_d_scores = expand_d_scores(d_scores,num_conditions,evt_list,subj_group)

   new_d_scores = [];
   num_in_grp = [0 num_conditions*subj_group];

   for grp_idx = 1:length(subj_group)
      num_subjs = subj_group(grp_idx);
      first = sum(num_in_grp(1:grp_idx)) + 1;
      last = sum(num_in_grp(1:(grp_idx+1)));
      tmp_evt_list = evt_list(first:last);
      tmp_d_scores = d_scores(((grp_idx-1)*num_conditions+1):grp_idx*num_conditions,:);
      tmp_d_scores = tmp_d_scores(tmp_evt_list,:);	% expand and reorder
      new_d_scores = [new_d_scores;tmp_d_scores];
   end

   return;