fmri_pls_analysis.m

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

            st_datamat = grp_datamat;
            saved_info = [saved_info, ' st_datamat'];
        end

     else

        if (num_boot > 0),
           num_rep = length(st_evt_list) / num_conditions;

%           boot_progress = rri_progress_ui('initialize');

           if isempty(subj_group)
              [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
                 = rri_boot_check(num_rep, 1, num_boot, 0, ...
                   for_batch);
%                   boot_progress, for_batch);
           else
              [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
                 = rri_boot_check(subj_group, num_conditions, num_boot, 0, ...
                   for_batch);
%                   boot_progress, for_batch);
           end

           num_boot = new_num_boot;
        else
           min_subj_per_group=[];is_boot_samples=[];boot_samples=[];new_num_boot=[];
        end;

        [brainlv,s,designlv,b_scores,d_scores,lvintercorrs,design, ...
            perm_result,boot_result,lv_evt_list] = ...
		fmri_taskpls_norotate(st_datamat,contrasts,num_conditions,...
			st_evt_list,num_boot,num_perm,subj_group, ...
         min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

        saved_info=['brainlv s designlv perm_result boot_result st_coords ', ...
  	         'st_dims lv_evt_list st_win_size st_voxel_size st_origin ', ...
                 'SessionProfiles ContrastFile b_scores d_scores design ', ...
		 'subj_group num_conditions cond_name cond_selection ', ...
		 'num_subj_lst subj_name lvintercorrs design ', ...
		 'session_files_timestamp datamat_files_timestamp create_ver'];

        if save_datamat & ~isempty(brainlv)
            first = 1;
            last = 0;
%            grp_datamat = [];
            for g = 1:length(num_subj_lst)
               last = last + num_conditions*num_subj_lst(g);
               [tmp idx] = sort(st_evt_list(first:last));
               tmp = st_datamat(first:last,:);
%               grp_datamat = [grp_datamat; tmp(idx,:)];
               grp_datamat{g} = tmp(idx,:);
               first = last + 1;
            end;
        
            st_datamat = grp_datamat;
            saved_info = [saved_info, ' st_datamat'];
        end

     end

  end;


   %  Either used "single" in analysis or had "single" in datamat
   %
   if singleanalysis | singledatamat
      singleprecision = 1;
   else
      singleprecision = 0;
   end
   saved_info = [saved_info, ' singleprecision'];


  %  Requested by Nancy, because they no longer use "dp" or "dprob" to evaluate
  %  designlv
  %
  if ~isempty(perm_result) & ~isbehav
     perm_result = rmfield(perm_result, 'dp');
     perm_result = rmfield(perm_result, 'designlv_prob');
  end

  if exist('progress_hdl','var') & ishandle(progress_hdl)

     progress_bar = getappdata(gcf,'ProgressBar');
     set(progress_bar,'Position',[0 0 1 1],'Visible','on');

  end

  if ~for_batch
     disp('RunPLS is done ...');
  end

  if isempty(brainlv)
     done = 1;
  else
     done = 0;
  end

  while ~done
    try
       eval(['save ''' resultFile ''' ' saved_info]);
       done = 1;
    catch
       if findstr('BfMRIsession.mat', fn)
          [result_file,result_path] = uiputfile('*BfMRIresult.mat', ...
			'Can not write file, please try again');
       else
          [result_file,result_path] = uiputfile('*fMRIresult.mat', ...
			'Can not write file, please try again');
       end

       if isequal(result_file,0)		% Cancel was clicked
          resultFile = [];
          msg1 = ['WARNING: No file is saved.'];
%          uiwait(msgbox(msg1,'Uncompleted','modal'));
          return;
       else
          resultFile = fullfile(result_path,result_file);
       end;
    end
  end

  return; 					% fmri_pls_analysis


%-------------------------------------------------------------------------
function contrasts = load_contrast_file(contrastFile)

  load(contrastFile);

  num_contrasts = length(pls_contrasts);
  num_conditions = length(pls_contrasts(1).value);

  contrasts = zeros(num_conditions,num_contrasts);
  for i=1:num_contrasts,
     contrasts(:,i) = pls_contrasts(i).value';
  end;

  return;					% get_contrast


%-------------------------------------------------------------------------
function  [SessionProfiles,ContrastFile,num_perm,single_subj] = options_query()

  %  get the session profiles
  %
  num_session = 0;
  not_done = 1;
  while (not_done)
     msg = [ '\nEnter the name of file contains the session information ', ...
             '(Press <Enter> when done): '];
     sessionFile = input(msg,'s');

     if (isempty(sessionFile))
        not_done = 0;
     else
        if sessionFile(1) ~= '/',
            sessionFile = sprintf('%s/%s',pwd,sessionFile);
        end;
        num_session = num_session + 1;
        SessionProfiles{num_session} = sessionFile;
     end;
  end;


  msg = ['\nEnter the contrast file: ', ...
         '\n   (Type "HELMERT" for Helmert matrix or ', ...
         '\n    press <Enter> to use deviation from grand mean) '];
  ContrastFile = input(msg,'s');

  %  get the number of iterations for the permutation test 
  %
  num_perm = input('\nNumber of permutations: ');


  %  single subject analysis ?
  %
  single_subj_ans=input('\nDoes the data come from single subject? [Y/N] ','s');
  if (upper(single_subj_ans) == 'Y')
     single_subj = 1;
  else
     single_subj = 0;
  end;

  return;					% options_query


%-------------------------------------------------------------------------
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;
     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_st_datamat, new_st_coords, st_dims, num_conditions, new_evt_list, ...
	    win_size, voxel_size, origin, subj_group, ...
			subj_name, curr_conditions, num_behav_subj, ...
			newdata_lst, num_subj_lst ] = ...
    			concat_st_datamat2(singleanalysis,SessionProfiles,progress_hdl,ContrastMethod, ...
			posthoc, behavdata, cond_selection)
%
%
  newdata_lst = {};
  num_subj_lst = [];

  new_st_datamat = [];			%  stacked datamat
  new_st_coords = [];
  st_dims = [];
  num_conditions = [];
  new_evt_list = [];
  new_evt_list_lst = {};
  win_size = [];
  voxel_size = [];
  origin = [];
  subj_group = [];
  subj_name = {};
  curr_conditions = {};
  num_behav_subj = [];
  num_behavdata_col = size(behavdata,2);

  msg = sprintf('Merging datamats ...');
  ShowProgress(progress_hdl,msg);

  num_groups = length(SessionProfiles);

  profile_list = [];
  session_group = zeros(1,num_groups);
  for i=1:num_groups,
     session_group(i) = length(SessionProfiles{i});
     profile_list = [profile_list; SessionProfiles{i}];
  end;
  num_profiles = length(profile_list);

  st_info = cell(1,num_profiles);

  %  get the coords ...
  %
  total_evts = 0;
  subj_group = zeros(1,num_groups);
  cnt = 0;

  fn = SessionProfiles{1}{1};

  for i=1:num_groups,
    for j=1:session_group(i),
       cnt = cnt+1;

       ShowProgress(progress_hdl,cnt/(num_profiles*10));
       msg = sprintf('Loading datamat: %d out of %d',cnt,num_profiles);
       ShowProgress(progress_hdl,msg);

       sessionFile = profile_list{cnt};
       load(sessionFile);
       curr_conditions = session_info.condition;

        %  get the st_datamat and st_evt_list 
        %
        datamat_prefix = session_info.datamat_prefix;

        if findstr('BfMRIsession.mat', fn)
           st_datamatFile = sprintf('%s_BfMRIdatamat.mat',datamat_prefix);
        else
           st_datamatFile = sprintf('%s_fMRIdatamat.mat',datamat_prefix);
        end

%        st_datamatFile = fullfile(session_info.pls_data_path,st_datamatFile);

        curr = pwd;
        if isempty(curr)
           curr = filesep;
        end

        st_datamatFile = fullfile(curr,st_datamatFile);
   
	try
           warning off;

           load(st_datamatFile,'st_coords','st_dims','st_evt_list', ...
			'st_win_size','st_voxel_size','st_origin');
           warning on;
	catch
	   disp(sprintf('ERROR: cannot open data file: %s',st_datamatFile));
	   new_st_datamat = [];
	   return;   
	end;

        st_info{cnt}.sessionFile = sessionFile;
        st_info{cnt}.datamatFile = st_datamatFile;
        st_info{cnt}.num_cond = session_info.num_conditions;
        st_info{cnt}.coords = st_coords;
        st_info{cnt}.dims = st_dims;
        st_info{cnt}.evt_list = st_evt_list;
        st_info{cnt}.win_size = st_win_size;
        st_info{cnt}.voxel_size = st_voxel_size;
        st_info{cnt}.origin = st_origin;

        num_evts = length(st_evt_list);
        total_evts = total_evts + num_evts;
        subj_group(i) = subj_group(i) + num_evts/length(curr_conditions);
   
        if (cnt > 1),	% make sure the st_datamat are compatible
           if (st_info{cnt-1}.dims ~= st_info{cnt}.dims),
              msg = 'The datamats have different volume dimension.';
              ShowProgress(progress_hdl,['ERROR: ', msg]);
disp(msg);
              return;
           end;
   
           if (st_info{cnt-1}.win_size ~= st_info{cnt}.win_size),
              msg = 'The datamats have different window size.';
              ShowProgress(progress_hdl,['ERROR: ', msg]);
disp(msg);
              return;
           end;

           if (st_info{cnt-1}.dims ~= st_info{cnt}.dims),
              msg = 'The datamats have different volume dimension.';
              ShowProgress(progress_hdl,['ERROR: ', msg]);
disp(msg);
              return;
           end;

           if (st_info{cnt-1}.voxel_size ~= st_info{cnt}.voxel_size),
              msg = 'The datamats have different voxel size.';
              ShowProgress(progress_hdl,['ERROR: ', msg]);
disp(msg);
              return;
           end;
   
           if ~isequal(curr_conditions,prev_conditions)
              msg='The datamats are created from different conditions.';
              ShowProgress(progress_hdl,['ERROR: ', msg]);
disp(msg);
              return;
           end;
        end;  % (cnt > 1)
   
        prev_conditions = curr_conditions;
   
    end;
  end;

  clear st_datamat;

  num_conditions = length(curr_conditions);  
  st_dims = st_info{1}.dims; 

  %  determine the common coords ...
  %
  m = zeros(1,prod(st_info{1}.dims));
  for i=1:num_profiles
     m(st_info{i}.coords) = m(st_info{i}.coords) + 1;
  end;
  new_st_coords = find(m == num_profiles);

  if isempty(new_st_coords)
     disp('ERROR: no common coords among datamats!');
     new_st_datamat = [];
     return;
  end

  %  stack the st_datamat together
  %
  win_size = st_info{1}.win_size;
  voxel_size = st_info{1}.voxel_size;

  if isempty(st_info{1}.origin) | isequal(st_info{1}.origin,[0 0 0]),
     origin = floor(st_dims([1 2 4])/2); 
  else
     origin = st_info{1}.origin;
  end;

  num_voxels = length(new_st_coords);
  new_st_datamat = zeros(total_evts,num_voxels*win_size);
  first_row = 1;

  first_cond_order = [];

  %  go through each subject, which is represented by each profile
  %
%  for i=1:num_profiles,

  tmp_new_st_datamat = [];
  cnt=0;

  for i=1:num_groups,

    grp_tmp_new_st_datamat = [];
    grp_tmp_new_evt_list = [];
    grp_tmp_new_evt_list = [];
    grp_first_cond_order = [];

    for j=1:session_group(i),
       cnt = cnt+1;

       ShowProgress(progress_hdl,cnt/(num_profiles*10)+1/10);

       load(st_info{cnt}.sessionFile);
       load(st_info{cnt}.datamatFile);
       Subj = session_info.datamat_prefix;

       if singleanalysis
          st_datamat = single(st_datamat);
       else
          st_datamat = double(st_datamat);
       end

       coord_idx = find( m(st_info{cnt}.coords) == num_profiles ); 
       nr = length(st_info{cnt}.evt_list);	% number of runs
       nc = length(st_info{cnt}.coords);
       last_row = nr + first_row - 1; 


if ContrastMethod == 4 | ContrastMethod == 5

       %  stack behavdata and get behavmask (re-order for each session file
       %  to make it 'each condition in each run (yes, reversed)'
       %
       behavmask = 1:size(st_datamat,1);

       %  nrr could be nr, depend on 'across run' or 'within run'
       %
       nrr = size(st_datamat,1) / num_conditions;

       behavmask = reshape(behavmask, [nrr, num_conditions]);
       behavmask = behavmask';
       behavmask = reshape(behavmask, [size(st_datamat,1),1]);

end


       this_subj_order = zeros(1, nr);
       first_cond = 1;
       jj = 1;

       %  get first_cond of each run
       %
       while first_cond <= nr
          this_subj_order(first_cond) = 1;
          first_cond = first_cond + num_conditions;

          if exist('SingleSubject','var') & SingleSubject == 1
             subj_name = [subj_name, {[Subj, '_Onset', num2str(jj)]}];
          else