function [Cond] = rd_taskcorr(varargin);
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      

% Correlations between task and global brain mean or motion.
% FORMAT Cond = rd_taskcorr;
%
% FORMAT Cond = rd_taskcorr('Global');
% Computes the correlations between the different task 
% conditions and the global brain mean intensity.
%
% FORMAT Cond = rd_taskcorr('Motion');
% Computes the correlations between the different task 
% conditions and the motion parameters estimated during
% the realignment stage.
%
% FORMAT Cond = rd_taskcorr('Global','Motion');
% FORMAT Cond = rd_taskcorr('Motion','Global');
% Computes the correlations with both the global brain
% mean and the motion parameters.
%
% Returns a table of the correlations and their corresponding
% significance values. The correlations are computed over
% the averaged timcourse during the peri-simulus period (or
% Cond.min_elapse). Please see the Cond structure description
% below for more information.
% 
% Cond    - structure containing information for each condition
% .name   - name of condition
% .iC     - corresponding column(s) in the design matrix xX.X
% .iS     - corresponding condition indeces in each Sess
% .iSC    - session specific indeces for this condition
% .X      - task variates from xX.X collapsed across .iC columns
% .g      - the average global brain during the peri-stimulus time
% .G      - the average task variate during the peri-stimulus time
% .rGg    - correlation (Pearson's R) between task and global mean
% .t      - t statistic corresponding to .rGg
% .p      - p value giving the significance of .rGg (two-tailed)
% .onsec  - the onsets for this condition (in sec)
% .onidx  - the onsets for this condition (in scans)
% .num_scans  - (a guess about the) duration of the PST in scans
% .min_elapse - the minimum time between consecutive presentation
%               of trials or blocks of this condition type (it is
%               used to make a guess about the duration of the PST)
% 
% If the 'Motion' option is selected, the following fields will
% be added to Cond:
%
% .RP     - a 6xPST matrix of the average for each of the 
%           realignment parameters during the peri-stimulus time
% .RPname - name of each realignment parameter and the averages
% .RPr    - correlations (Pearson's R) between task and the 6 
%           realignment parameters, the average of translations, 
%           the  average of rotations, and the average of all 6
%           realignment parameters
% .RPt    - t statistics corresponding to .RPr
% .RPp    - p values giving the significance of .RPt (two-tailed)
%
% ---------------------------------------------------------------
% Compatible with SPM99 files. 
% Uses SPMcfg.mat to load parameters.
% 
% Condition names across different sessions are matched ignoring 
% letter case. The order in which conditions have been entered in 
% each session is not important.
% ---------------------------------------------------------------
% @(#) rd_taskcorr.m 1.1      Kalina Christoff         2000-06-10 

Option=[];

% if no input, ask
% -----------------------------------
if nargin == 0;
     fprintf('\n\nPlease select options from below.\n\n');  

     fprintf('Compute correlations with global brain mean? Type [y] or ');
     if strcmpi(input('enter to skip :','s'),'y'); Option={'Global'}; end
                
     fprintf('Compute correlations with motion parameters? Type [y] or');
     if strcmpi(input('enter to skip :','s'),'y'); 
        Option={char(Option),'Motion'};
     end
end


% if input is given, assign
% -----------------------------------

if nargin ==1 | nargin ==2;
  if  find(strcmpi('Global',cellstr(varargin)))  ; Option={'Global'}; end
  if  find(strcmpi('Motion',cellstr(varargin))) ; Option={char(Option), 'Motion'}; end
end


% if more than 2 input argument, complain that there's too much information
% -------------------------------------------------------------------------
if nargin > 2; error('Too many input arguments. Please try again.'); end


% load things
% ==========================================================
SPMcfg_path = spm_get(1,'SPMcfg.mat',...
                     'Please select design specification file.');
load(SPMcfg_path);

% get # sessions
n_sess         = length(Sess);


% if the the option includes 'Motion', load the realignment params
% ================================================================
if find(strcmp(Option,'Motion'))

    RP_path = spm_get(n_sess,'.txt',...
		    'Select a motion param file for each session.');
    RP=[];

    for s=1:n_sess
       RP=[RP; spm_load(deblank(RP_path(s,:)))];
    end

    % compute the index for which params correspond to existing V-files
    % This is to deal with cases when not all V-files are entered in 
    % the statistical model
    % ==================================================================
    if length(VY) < size(RP,1)
    
	fprintf('\n\nIt appears that some V.img files were excluded from the ');
	fprintf('analysis.\nPlease wait while excluding them from the motion ');
        fprintf('parameters file.\n\n');
    
        % get only the indeces for the first session only
        % -----------------------------------------------
        excluded=str2num(VY(Sess{1}.row(1)).fname(end-6:end-4))-1;
        RPidx=[zeros(excluded,1); ones(length(Sess{1}.row),1)];
            fprintf('Assuming that the first  %d files ',excluded);
            fprintf('from session  1 were excluded.\n');
      
        % get the other sessions
        % -----------------------------------------------
            add=excluded;
        for s=2:n_sess
            excluded=str2num(VY(Sess{s}.row(1)).fname(end-6:end-4))-1;
            add=add+excluded;
            RPidx=[RPidx; zeros(excluded,1); ones(length(Sess{s}.row),1)];
                fprintf('Assuming that the first  %d files from session  %d ',...
                         excluded,s);fprintf('were excluded.\n');
        end;  clear excluded add;

	RP=RP(find(RPidx),:);
    
     % finish compensating for excluded V files
     % ========================================    
     end    

     if length(VY) ~= size(RP,1)
	    fprintf('\n\nSorry, the motion parameters still do not match the');
            fprintf('\nnumber of V.img files that were in the analysis.\n');
            fprintf('\nPlease, try again.\n\n');
	    break
     end

% finish loading and processing realignment params    
% ========================================================================
end
    
    
% find the correct condition index corresponding to the same name
% in each session, so that even conditions entered in different order
% for the different sessions would be averaged correctly
% ---------------------------------------------------------------------

% get the names of all conditions in the entire experiment
% ----------------
  cond_names = [];
  xX.iC         = [];    % the indeces in xX.X for hrf only

  for s=1:n_sess
    cond_names = [cond_names, Sess{s}.name];
    tmp=[];
    for i=1:length(Sess{s}.ind)
       tmp    = [tmp, Sess{s}.ind{i}(1)]; 
    end
    xX.iC = [xX.iC, Sess{s}.col(tmp)];
    clear tmp;
  end
  
% guess number of conditions
% --------------------------

name_idx = 1:length(xX.iC);

c=1; % start from condition one

while ~isempty(find(name_idx))
 
   Cond(c).name = cond_names(min(find(name_idx)));

	   % match the condition names, ignoring case
	   % ----------------------------------------
	   Cond(c).iC   = xX.iC(find(strcmpi(Cond(c).name,cond_names)))';
		
		  %-----------------------------------------------------------
	     if length(Cond(c).iC) < n_sess  % if this condition does not
                                        % appear in every session
		       s=1;
				 while s <= length(Cond(c).iC)
			
					i = s;
					while  Cond(c).iC(i) > max(Sess{i}.col)		
								Cond(c).iC(i+1:end+1,1)...
							   =Cond(c).iC(i  :end  ,1);
				  
	 						   Cond(c).iC(i)=0;
					  i = i+1;
					end
				   s = s+1;

				 end
			 
 		       Cond(c).iC(end+1:n_sess,1)=zeros;

		  end
		  % --------------------------------------------------------------

		  

			  
   Cond(c).iS  = zeros(n_sess,1);
   Cond(c).iSC = zeros(n_sess,1);
   cnt = 0; 

   % update name_idx with zeros where used already
   % ---------------------------------------------
   for s = 1:n_sess

    idx = find(xX.iC==Cond(c).iC(s));
    
    if ~isempty(idx)
      % Cond.iS is the indeces for this condition in each sessions
      % ----------------------------------------------------------
      Cond(c).iS(s) = name_idx(idx);

      % Zero out the condition indeces that have been used
      % --------------------------------------------------
      name_idx(idx) = 0;    

      % Cond.iSC is the session specific condition index
      % ------------------------------------------------
      Cond(c).iSC(s) = Cond(c).iS(s) - cnt;
    end
    
    cnt = cnt+length(Sess{s}.name);

   end
   
   c=c+1;                                   % update to next condition
end

clear cond_names s c name_idx idx cnt;

% compose the G matrix by collapsing xX.X across session
% ------------------------------------------------------
n_cond=length(Cond);

for c = 1:n_cond
      Cond(c).X = sum(xX.X(:,Cond(c).iC(find(Cond(c).iC))),2);  
end


% START LOOPING THROUGH CONDITIONS
% ================================
for c=1:n_cond


    % convert all Sess{i}.ons into secs from the beginning of the experiment
    % ----------------------------------------------------------------------
    pst_ons=[];        % the onsets in sec

    for s = 1:n_sess
       if Cond(c).iSC(s)~=0          % skip this session if empty
         pst_ons = [pst_ons; ...
                   Sess{s}.ons{Cond(c).iSC(s)} + ...
                   + (min(Sess{s}.row)-1)*xX.RT ...
                 ];
       end
    end

    % get the onsets and the offsets for the entire experiment, in sec
    % ----------------------------------------------------------------
    pst_ons_idx  = floor(pst_ons);
    pst_offs_idx = [(pst_ons_idx(2:end)-1);max(Sess{n_sess}.row)];

    % make a guess how long was the trial duration, i.e. 
    % take the minimum time between consecutive presentation of events 
    % from the same condition
    % ----------------------------------------------------------
    elapses    = pst_offs_idx-pst_ons_idx+1;          % in seconds


    % If there is only one repetition, plot timecourse till end 
    % of session. Otherwise, get the minimum of all lapses, except 
    % the last one (which is defined by where the last session ends).
    % --------------------------------------------------------------------
    if      length(elapses)==1, min_elapse = elapses;
    elseif  length(elapses) >1,
                                min_elapse = ...
                                min(elapses(1:end-1));   
                                % min time in seconds
    end

    num_scans  = round(min_elapse/xX.RT); % min time in scans

    % convert the pst_ons_idx  into scans, to use as indeces
    % ----------------------------------------------------------------------
    pst_ons_idx  = floor(pst_ons_idx/xX.RT)+1;        % indeces in scans
     
    % average across PST for
    % - the global brain variate g
    % - the task variate Cond(i).X
    % - the six realignment parameters
    % ----------------------------------------------------------------

	 for i=1:num_scans
    	     % for the last trial in the experiment, average only
             % over the available scans
             % ------------------------
	     if max(pst_ons_idx)+i-1 <= max(Sess{n_sess}.row);
     
               Cond(c).g(i)      =mean(xGX.rg   (pst_ons_idx+i-1));
	       Cond(c).G(i)      =mean(Cond(c).X(pst_ons_idx+i-1));

               if exist('RP')
                  Cond(c).RP(:,i)=mean(RP(pst_ons_idx+i-1,:),1);
               end

             elseif max(pst_ons_idx)+i-1 > max(Sess{n_sess}.row);

               Cond(c).g(i)=mean(xGX.rg   (pst_ons_idx(1:end-1)+i-1));
	       Cond(c).G(i)=mean(Cond(c).X(pst_ons_idx(1:end-1)+i-1));    

               if exist('RP')
                  Cond(c).RP(:,i)=mean(RP(pst_ons_idx(1:end-1)+i-1,1),1);
               end
    
             end    

         end


    % calculate correlations with task
    % =========================================================================
    Cond(c).rGg  = (corrcoef(Cond(c).g,Cond(c).G)*[1;0])'*[0;1]; 

    % test the significance of correlation (assumes normal distribution though)
    % -------------------------------------------------------------------------
    Cond(c).t    = (Cond(c).rGg*sqrt(num_scans-2))/...
                         sqrt(1-Cond(c).rGg^2);

    % t follows the t(n-2) distribution; compute 2-tailed significance
    % ------------------------------------------------------------
    if     Cond(c).t >= 0,
                Cond(c).p    = (1-tcdf(Cond(c).t,num_scans-2))*2;
    else        Cond(c).p    = tcdf(Cond(c).t,num_scans-2)*2;
    end



    % if the option includes 'Motion', compute correlations
    % =========================================================================
    if exist('RP')

            Cond(c).RPname = {'x_mm','y_mm','z_mm','pitch','roll','yaw',...
                              'transl','rotat','all'};

            Cond(c).RPr    = corrcoef([Cond(c).G;Cond(c).RP;...
                                       mean(Cond(c).RP(1:3,:));...
                                       mean(Cond(c).RP(4:6,:));...
     			               mean(Cond(c).RP(1:6,:));...
                                      ]');
                             Cond(c).RPr=Cond(c).RPr(1,2:10);

        % test the significance of correlation (assumes normal distribution though)
        % -------------------------------------------------------------------------
        Cond(c).RPt    = (Cond(c).RPr*sqrt(num_scans-2))./...
                            sqrt(1-Cond(c).RPr.^2);

	% t follows the t(n-2) distribution; compute 2-tailed significance
	% ------------------------------------------------------------
        for i=1:9
	    if     Cond(c).RPt(i) >= 0,
		     Cond(c).RPp(i) = (1-tcdf(Cond(c).RPt(i),num_scans-2))*2;
	    else     Cond(c).RPp(i)    = tcdf(Cond(c).RPt(i),num_scans-2)*2;
	    end
        end
     end
     %====================


    Cond(c).onsec      = pst_ons;       % onsets in sec
    Cond(c).onidx      = pst_ons_idx;   % onsets in scans
    Cond(c).num_scans  = num_scans;
    Cond(c).min_elapse = min_elapse;


    
	 
% END LOOPING THROUGH CONDITIONS
% ==============================
end

    % Print table with task correlations, if the option includes 'Global'
    % -----------------------------------------------------------------
    if find(strcmp(Option,'Global'))

	fprintf('\n\n---------------------------------------------------------------\n');
	fprintf('Correlations between global brain intensity and different task\n');
	fprintf('conditions, with corresponding significance values (two-tailed).\n');
	fprintf('---------------------------------------------------------------\n\n');   

	fprintf('%20s%17s%15s\n','Condition','Pearson''s R','p value');     

	for c=1:length(Cond)
	 fprintf('%20s%15.2f%15.2f\n',char(Cond(c).name), Cond(c).rGg, Cond(c).p);
	end

	fprintf('\n---------------------------------------------------------------\n');
	fprintf('                                                 @(#) RD 1.1 KC\n\n'); 
    end


    % Print table with motion correlations, if the option include 'Motion'
    % --------------------------------------------------------------------
    if find(strcmp(Option,'Motion')) & exist('RP')
        
        more on;
        fprintf('\n\n\n---------------------------------------------------------------------------------\n');
	fprintf('Correlations between motion parameters and different task ');
	fprintf('conditions.\nPearson''s R values are given together with corresponding ');
        fprintf('significance values\n(in parentheses, two-tailed).\n\n');
        fprintf('                       ''transl'' is the mean of x_mm, y_mm, and z_mm\n');
        fprintf('                        ''rotat'' is the mean of pitch, roll, and yaw\n');
        fprintf('                          ''all'' is the mean of all six parameters\n');
	fprintf('---------------------------------------------------------------------------------\n\n');   

        fprintf('%15s%8s%7s%7s%8s%7s%7s%8s%7s%7s\n\n','Condition',...
                'x_mm','y_mm','z_mm','pitch','roll','yaw','transl','rotat','all');

	for c=1:length(Cond)
	 fprintf('%15s%8.2f%7.2f%7.2f%8.2f%7.2f%7.2f%8.2f%7.2f%7.2f\n',...
                 char(Cond(c).name),Cond(c).RPr);
         fprintf('\t\t  (%1.2f) (%1.2f) (%1.2f)  (%1.2f) (%1.2f) (%1.2f)  (%1.2f) (%1.2f) (%1.2f)\n',Cond(c).RPp);
    	end

	fprintf('\n----------------------------------------------------------------------------------\n');
	fprintf('                                                                   @(#) RD 1.0 KC\n\n'); 
    end












