pet_analysis.m

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

%PET_ANALYSIS Apply PLS on the PET data based on the information saved in
%    the datamat file.
%
%    Usage: [resultFile, elapsed_time] = ...
%	pet_analysis(isbehav, datamat_files, num_perm, ...
%	num_boot, Clim)
%
%    see also PLS_FMRI_ANALYSIS
%

%   Called by pet_analysis_ui
%
%  INPUT:
%    isbehav - 1 if run Behavior PLS; 0 for Task PLS.
%    datamat_files - a cell array, one element per group.  Each element
%               in the array is another cell array contains the names of
%               session profiles for the group.
%    num_perm - number of permutations to be performed.
%    num_boot - number of bootstrap resampling to be performed.
%    Clim - upper limit of confidence interval estimated
%
%   OUTPUT FILE:
%         - file stores the information of the PLS result.
%
%   Created July 2001 by Wilkin Chau, Rotman Research Institute
%   Modified on 02-OCT-2002 by Jimmy Shen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [resultFile, elapsed_time] = pet_analysis(varargin)

  singledatamat = 0;		% init singledatamat to false

  isbehav = varargin{1};
  datamat_files = varargin{2};
  num_perm = varargin{3};
  num_boot = varargin{4};
  Clim = varargin{5};
  posthoc = varargin{6};
  cond_selection = varargin{7};
  behavname = varargin{8};
  behavdata = varargin{9};
  behavdata_lst = varargin{10};
  ContrastFile = varargin{11};
  iscontrast = varargin{12};
  ismean = varargin{13};
  save_datamat = varargin{14};
  ismultiblock = varargin{15};
  bscan = varargin{16};

  if (nargin > 16)
     output_file = varargin{17};
     for_batch = 1;
  else
     for_batch = 0;
  end;

  datamat_files_timestamp = datamat_files;

  for i = 1:length(datamat_files)
     tmp = dir(datamat_files{i});
     datamat_files_timestamp{i} = tmp.date;

        warning off;
        load(datamat_files{i}, 'singleprecision');
        warning on;

        if exist('singleprecision','var') & singleprecision
           singledatamat = 1;
        end
  end

  if exist('output_file','var') & ~isempty(output_file)
    resultFile = output_file;
  else

    load(datamat_files{1},'session_info');
    datamat_prefix = session_info.datamat_prefix;

    [result_file,result_path] = ...
	uiputfile([datamat_prefix '_PETresult.mat'],'Saving PLS Result');

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


   v7 = version;
   if str2num(v7(1))<7
      singleanalysis = 0;
   else
      singleanalysis = 1;
   end

   pc = computer;
   if singleanalysis & ( strcmp(pc,'GLNXA64') | strcmp(pc,'GLNXI64') | strcmp(pc,'PCWIN64') )
      quest = questdlg({'We detected that you are running MATLAB on a 64-bit system. According to MATLAB Bug Report ID 268001, we have to convert data to double precision for Intel based system.' '' 'Is this Intel 64-bit machine?' ''}, 'Choose','No','Yes','Don''t know','Don''t know');

      if ~strcmp(quest,'No')
         singleanalysis = 0;
      end
   end;


  progress_hdl = rri_progress_ui('initialize');

  % load the session info and datamat
  %
  [behavdata_lst, newdata_lst, newcoords, dims, num_cond_lst, ...
	num_subj_lst, subj_name_lst, voxel_size, origin, ...
	behavname, behavdata] = pet_get_common(datamat_files,  ...
	cond_selection, behavname, behavdata, ...
	behavdata_lst, progress_hdl);

   if isempty(newcoords)
      disp('ERROR: There are no common voxels found.');
      return;
   else
      for i = 1:length(newdata_lst)
         if isempty(newdata_lst{i})
            disp('ERROR: Merged datamat is empty.');
            return;
         end
      end
   end


   for grp=1:length(newdata_lst)
      if singleanalysis
         newdata_lst{grp} = single(newdata_lst{grp});
      else
         newdata_lst{grp} = double(newdata_lst{grp});
      end
   end


  %  start PLS Run...
  %

  perm_result = [];
  boot_result = [];
  create_ver = plsgui_vernum;

  if(ismultiblock)

     isbehav = 2;

     rri_progress_ui(progress_hdl, 'Running Multiblock PLS', 'Running Multiblock PLS ...');

     ibehavdata_lst = behavdata_lst;

     if (num_boot > 0)
%        boot_progress = rri_progress_ui('initialize');

        [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
           = rri_boot_check(num_subj_lst, num_cond_lst(1), num_boot, 0, ...
             for_batch);
%             boot_progress, for_batch);

        num_boot = new_num_boot;
     end

     if (num_perm > 0) | (num_boot == 0)
        [brainlv,s,designlv,behavlv,brainscores,designscores,behavscores, ...
		lvcorrs, origpost, perm_result, datamatcorrs_lst, ...
		b_scores,behav_row_idx,behavdata_lst] = ...
		pet_multiblock_perm(ibehavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_perm,isbehav,posthoc,bscan);
     end

     if (num_boot > 0)
        if num_perm == 0, origpost = []; end;
        [brainlv2,s2,designlv2,behavlv2,brainscores2,designscores2, ...
		behavscores2,lvcorrs2, boot_result, datamatcorrs_lst2, ...
		b_scores2,behav_row_idx2,behavdata_lst2] = ...
		pet_multiblock_boot(ibehavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_boot,isbehav,Clim, ...
		min_subj_per_group,is_boot_samples,boot_samples, ...
		new_num_boot,bscan);

        if num_perm == 0
           brainlv = brainlv2;
           s = s2;
           designlv = designlv2;
           behavlv = behavlv2;
           brainscores = brainscores2;
           designscores = designscores2;
           behavscores = behavscores2;
           lvcorrs = lvcorrs2;
           perm_result = [];
           datamatcorrs_lst = datamatcorrs_lst2;
           b_scores = b_scores2;
	   behav_row_idx = behav_row_idx2;
           behavdata_lst = behavdata_lst2;
        end
     end

     saved_info=['''brainlv'', ''s'', ''designlv'', ''behavlv'', ''brainscores'', ', ...
		'''designscores'', ''behavscores'', ''lvcorrs'', ''origpost'', ',...
		'''perm_result'', ''boot_result'', ''datamatcorrs_lst'', ', ...
		'''newcoords'', ''cond_selection'', ''dims'', ''b_scores'', ', ...
		'''voxel_size'', ''origin'', ''bscan'', ''behavname'', ', ...
		'''num_cond_lst'', ''num_subj_lst'', ''subj_name_lst'', ', ...
		'''behavdata_lst'', ''ismultiblock'', ''datamat_files'', ', ...
		'''datamat_files_timestamp'', ''create_ver'''];

     if save_datamat
         saved_info = [saved_info, ', ''newdata_lst'''];
     end

  elseif(isbehav)

     rri_progress_ui(progress_hdl, 'Running Behavior PLS', 'Running Behavior PLS ...');

     if (num_boot > 0)
%        boot_progress = rri_progress_ui('initialize');

        [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
           = rri_boot_check(num_subj_lst, num_cond_lst(1), num_boot, 0, ...
             for_batch);
%             boot_progress, for_batch);

        num_boot = new_num_boot;
     end

     if (num_perm > 0) | (num_boot == 0)
        [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
		origpost, perm_result, datamatcorrs_lst] = ...
		pet_analysis_perm(behavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_perm,isbehav,posthoc);
     end

     if (num_boot > 0)
        if num_perm == 0, origpost = []; end;
        [brainlv2,s2,behavlv2,brainscores2,behavscores2,lvcorrs2, ...
		boot_result, datamatcorrs_lst2] = ...
		pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_boot,isbehav,Clim, ...
		min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

        if num_perm == 0
           brainlv = brainlv2;
           s = s2;
           behavlv = behavlv2;
           brainscores = brainscores2;
           behavscores = behavscores2;
           lvcorrs = lvcorrs2;
           perm_result = [];
           datamatcorrs_lst = datamatcorrs_lst2;
        end
     end

     saved_info=['''brainlv'', ''s'', ''behavlv'', ''brainscores'', ', ...
		'''behavscores'', ''lvcorrs'', ''origpost'', ',...
		'''perm_result'', ''boot_result'', ''datamatcorrs_lst'', ', ...
		'''newcoords'', ''cond_selection'', ''dims'', ', ...
		'''voxel_size'', ''origin'', ''behavname'', ', ...
		'''num_cond_lst'', ''num_subj_lst'', ''subj_name_lst'', ', ...
		'''behavdata_lst'', ''datamat_files'', ', ...
		'''datamat_files_timestamp'', ''create_ver'''];

     if save_datamat
         saved_info = [saved_info, ', ''newdata_lst'''];
     end

  elseif(iscontrast)			% contrast analysis

     rri_progress_ui(progress_hdl, 'Running No-Rotate PLS', 'Running No-Rotate PLS ...');

     if isnumeric(ContrastFile)
        design = ContrastFile;
     else
        design = load(ContrastFile);
     end

     if (num_boot > 0)
%        boot_progress = rri_progress_ui('initialize');

        [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
           = rri_boot_check(num_subj_lst, num_cond_lst(1), num_boot, 0, ...
             for_batch);
%             boot_progress, for_batch);

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

     [brainlv,s,designlv,brainscores,designscores,lvintercorrs,design, ...
	perm_result,boot_result] = rri_taskpls_norotate(newdata_lst,design, ...
	num_subj_lst,num_cond_lst,num_boot,num_perm, ...
	min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

     saved_info=['''brainlv'', ''s'', ''designlv'', ''brainscores'', ', ...
		'''designscores'', ''lvintercorrs'', ''design'', ', ...
		'''perm_result'', ''boot_result'', ', ...
		'''newcoords'', ''cond_selection'', ''dims'', ', ...
		'''voxel_size'', ''origin'', ', ...
		'''num_cond_lst'', ''num_subj_lst'', ''subj_name_lst'', ', ...
		'''datamat_files'', ''datamat_files_timestamp'', ', ...
		'''create_ver'''];

     if save_datamat
         saved_info = [saved_info, ', ''newdata_lst'''];
     end

  else					% deviation analysis

     rri_progress_ui(progress_hdl, 'Running Task PLS', 'Running Task PLS ...');

     if (num_boot > 0)
%        boot_progress = rri_progress_ui('initialize');

        [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
           = rri_boot_check(num_subj_lst, num_cond_lst(1), num_boot, 1, ...
             for_batch);
%             boot_progress, for_batch);

        num_boot = new_num_boot;
     end

     if (num_perm > 0) | (num_boot == 0)
        [brainlv,s,designlv,brainscores,designscores,lvcorrs, ...
		origpost, perm_result, datamatcorrs_lst] = ...
		pet_analysis_perm(behavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_perm,isbehav,posthoc);
     end

     if (num_boot > 0)
        [brainlv2,s2,designlv2,brainscores2,designscores2,lvcorrs2, ...
		boot_result, datamatcorrs_lst2] = ...
		pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,...
		num_subj_lst,num_boot,isbehav,Clim, ...
		min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

        if num_perm == 0
           brainlv = brainlv2;
           s = s2;
           designlv = designlv2;
           brainscores = brainscores2;
           designscores = designscores2;
           lvcorrs = lvcorrs2;
           perm_result = [];
           datamatcorrs_lst = datamatcorrs_lst2;
        end
     end

     saved_info=['''brainlv'', ''s'', ''designlv'', ''brainscores'', ', ...
		'''designscores'', ', ...
		'''perm_result'', ''boot_result'', ', ...
		'''newcoords'', ''cond_selection'', ''dims'', ', ...
		'''voxel_size'', ''origin'', ', ...
		'''num_cond_lst'', ''num_subj_lst'', ''subj_name_lst'', ', ...
		'''datamat_files'', ''datamat_files_timestamp'', ', ...
		'''create_ver'''];


     if save_datamat
         saved_info = [saved_info, ', ''newdata_lst'''];
     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'''];


  %  save results
  %
  msg = 'Saving to the disk ...';

  if exist('progress_hdl','var') & ishandle(progress_hdl)
     rri_progress_ui(progress_hdl, 'Save', msg);
  end

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

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

  while ~done
    try
       eval(['save(''', resultFile, ''',' saved_info,');']);
       done = 1;
    catch
       [result_file,result_path] = uiputfile('*PETresult.mat', ...
			'Can not write file, please try again');

       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; 					% pet_analysis