function roi_deconvolve_process(Design, IRF_LENGTH, CON_INT)
% --------------------------------------------------------
% Main processing loop



% Defaults 
% ---------------
tol = 1.0e-10;      % effective zero for machine accuracy
debug = 0;

% Design is a vector of structs, one for each running of the
% roi_deconvolve code, i.e., one element for each subject X ROI
% combination.
% Each element has the following fields:
%   - roi_fname: full-path filename for ROI file
%   - roi_name: extensionless ROI name, cleaned up and ready
%   for saving/labeling purposes
%   - XYZvox: the coordinates, in voxels, of this ROI (3 x
%   n)
%   - XYZmm: the coordinates, in mm, of this ROI (3 x n)
%   - subj_name: subject name, cleaned and suitable for
%   labeling/saving
%   - resdir: full-path results directory for this subject
%   - prepend: string to be prepended to front of output
%   filenames; can be empty.
%   - outdir: full-path directory in which all output files
%   are placed
%   - des_fname: full-path filename for the appropriate
%   design file - SPMcfg.mat for SPM99, SPM.mat for SPM2
%   - SPMver: version of SPM used for this subject.  Either
%   2 or 99 - integer.
%   - nsess: number of sessions for the experiment
%   - TR: TR for the experiment
%   - n_lags: round(IRF_LENGTH / TR) + 1 = the number of
%   timepoints in each condition-specific impulse function 
%   that will be estimated.
%   - omitted: cell array (nsess x 1 long) of session-specific
%   scan indexes to be omitted from analysis.
%   - V: vector of image files structures, each with
%   following fields:
%       - iS: session index for this V
%       - G: global mean for this V
%       - vol: memory-mapped vol struct (see spm_vol for
%       details)
%   - Cond: vector of condition structures, each with
%   following fields:
%       - name: name of condition
%       - basemodel_flag: 1 if condition is to be included
%       in baseline model and hence have its effects removed
%       from availability for independent testing; 0 otherwise.
%       Currently disabled - always 0.
%       - ons: vector of onsets (in TRs)
%       - dur: vector of durations (in TRs)
%       - stim: n(images) x 1 vector containing 1 at every
%       timepoint this stimulus is on and 0 everywhere else.
%   - base_flags: a structure itself, which contains a variety of
%   logical flags (and a couple integer values), signaling 
%   different aspects of the chosen baseline design:
%       - base_timeflag: 0 = steady baseline; 1 =
%       time-varying
%       - base_condflag: 0 = baseline is mean of timecourse;
%       1 = baseline is mean of condition
%       - base_cond_idx: 0 if baseline is mean of condition;
%       otherwise integer specifying condition number
%       - base_region: 1 = same ROI is baseline region;
%       2 = grand mean of volume is baseline; 3 = mean of
%       user-specified mask is baseline; 4 = mean of
%       user-specified inverse mask is baseline
%       - base_mask: empty unless base_region_idx = 3; then
%       full-path filename of specified mask
%       - base_neg: empty unless base_region_idx = 4; then
%       full-path filename of negative mask (voxels to
%       exclude)
%       - base_template: empty unless base_region_idx = 4;
%       then full-path filename of mask template (voxels
%       from which inverse mask are excluded)
%   - preproc_flags: a structure itself, which contains a
%   variety of logical flags and a couple other things,
%   signaling different aspects of the chosen temporal model of the
%   data:
%       - trim: 0 = no trim; 1 = trim
%       - trim_amt: threshold of standard deviations beyond
%       mean of session above or below which values will be
%       trimmed.
%       - detr: 0 = no detrending; 1 = linear detrending
%       - lowpass: 0 = no low-pass filtering; 1 = low-pass
%       filtering
%       - low_choice: 'Gaussian' or 'hrf'; empty if lowpass
%       = 0
%       - low_param: 0 unless low_choice is 'Gaussian';
%       then integer specifying FWHM in secs of low-pass
%       filter kernel
%       - highpass: 0 = no high-pass filtering; 1 =
%       high-pass filtering
%       - high_choice: 'none' if highpass is 0, otherwise
%       'specify'
%       - high_param: 0 if highpass is 0; otherwise, integer
%       specifying period of interest in secs; filter
%       attenuates effects with periods longer than this
%       number.
%   - Contrasts: vector of structures specifying what linear
%   contrasts, if any, are to be performed.  May be empty.
%   Has the following fields if not:
%       - name: label for contrast
%       - c: value of contrast, set up as array of n x p
%       size, where p is the number of parameters and n is
%       the number of linear constraints used.  Each row is
%       one linear constraint.


CON_COEFF = norminv([((1 - (CON_INT / 100)) / 2) (1 - ((1 - (CON_INT / 100)) / 2))],0,1);
CON_COEFF = CON_COEFF(2);
                    % confidence coefficient for desired
                    % confidence interval on IRFS
                    % 1.96 corresponds to 95%; 1.645 would
                    % correspond to 90%.
                    % norminv returns positive and negative
                    % values; we only need positive.
ROI_vector = [];


% Begin loop over subject x ROIs
% ------------------------------------------
for i = 1:length(Design)
    
CurDes = Design(i);

% ROI is a struct containing all the relevant information
% about the current analysis in question, with the following
% fields:
%   - name: combined prepend, subj_name and roi_name from
%   Des - cleaned for saving.
%   - XYZvox: the coordinates, in voxels, of this ROI (3 x
%   n)
%   - XYZmm: the coordinates, in mm, of this ROI (3 x n)
%   - baseline: either single value or timecourse of values;
%   the baseline value which the signal timecourse is
%   divided by to get our percent signal timecourse
%   - rawY: the raw signal timecourse, before conversion to
%   percent signal change - n x 1 vector, where n is number
%   of images
%   - rawY_perc: the raw percent signal timecourse, before
%   filtering/preprocessing - n x 1 vector
%   - Y: the actual percent signal timecourse that's the
%   input into our deconvolution model - n x 1 vector
%   - X: the full X matrix for our analysis - block effects,
%   any detrending effects desired, all time-lagged versions
%   of all condititions.
%   - base_X: the baseline model X matrix - same as X, but
%   containing only block and linear effects, as well as any
%   conditions specified as baseline conditions.
%   - full_xtxinv: (X'X)^-1 for the full model
%   - full_xtxinvxt: ((X'X)^-1)X' for the full model
%   - base_xtxinv: (X'X)^-1 for the baseline model - block
%   effects, linear detrend if asked for, any condition
%   effects that have been specified
%   - base_xtxinvxt: ((X'X)^-1)X' for the baseline model
%   - Cond: vector of structures co-indexed with the one in
%   Design; this one specifies the output by condition and
%   has the following fields:
%       - name: condition name (cleaned for labeling)
%       - basemodel_flag: 1 if condition is to be included
%       in baseline model and hence have its effects removed
%       from availability for independt testing; 0 otherwise.
%       - irf: impulse response function for this condition.  By
%       default, is 32 seconds long, with one element per
%       TR.  Can change length in default IRF_LENGTH above.
%       - redcd_X: the X matrix for the reduced model for
%       this condition, i.e., the X matrix with all columns for
%       this condition removed (empty if basemodel_flag==1)
%       - redcd_xtxinv: (X'X)^-1 for the reduced model (X
%       matrix without this condition) (empty if
%       basemodel_flag==1)
%       - redcd_xtxinvxt: ((X'X)^-1)X' for the reduced model (empty if
%       basemodel_flag==1)
%   - Contrasts: vector of structures specifying what linear
%   contrasts, if any, are to be performed.  May be empty.
%   Has the following fields if not:
%       - name: label for contrast
%       - c: value of contrast, set up as array of n x p
%       size, where p is the number of parameters and n is
%       the number of linear constraints used.  Each row is
%       one linear constraint.
%       - cxtxinvct: C((X'X)^-1)C' for this contrast.
%       - constant: constant matrix used for further
%       analysis - see below for details of construction.

% initialize ROI struct for this iteration
ROI = struct('name', [CurDes.prepend CurDes.subj_name '_' CurDes.roi_name], ...
    'XYZvox', CurDes.XYZvox, 'XYZmm', CurDes.XYZmm, 'outdir', CurDes.outdir, 'baseline', [], ...
    'rawY', zeros(length(CurDes.V.vol),size(CurDes.XYZmm, 2)), 'Cond', [], 'Contrasts', []);

ROI.name = spm_str_manip(ROI.name, 'ud');
ROI.roi_name = CurDes.roi_name;
ROI.subj_name = CurDes.subj_name;

vinv_data = inv(CurDes.V.vol(1).mat);  %matrix converting mm->voxels
ROI.XYZvox = vinv_data(1:3,:)*[ROI.XYZmm; ones(1,size(ROI.XYZmm,2))];
ROI.XYZvox = round(ROI.XYZvox);
            
for j = 1:length(CurDes.Cond)
    ROI.Cond(j).name = CurDes.Cond(j).name;
    ROI.Cond(j).name = spm_str_manip(ROI.Cond(j).name, 'ud');
    %spm_str_manip('u') does string cleaning - leaves it
    %with no spaces, minuses, invalid chars, etc.
    ROI.Cond(j).basemodel_flag = CurDes.Cond(j).basemodel_flag;
    ROI.Cond(j).sess = CurDes.Cond(j).sess;
end

for j = 1:length(CurDes.Contrasts)
    ROI.Contrasts(j) = CurDes.Contrasts(j);
    ROI.Contrasts(j).name = spm_str_manip(ROI.Contrasts(j).name, 'ud');
end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2.
% Extract raw timecourse from desired region and get baseline.
% -----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%ROI.Y = [
%    100.46;...
%    103.14;...
%    112.46;...
%    114.68;...
%    118.93;...
%    108.30;...
%    109.71;...
%    117.30;...
%    119.24;...
%    117.04;...
%    117.06;...
%    118.47;...
%    126.47;...
%    118.81;...
%    120.54;...
%    113.44;...
%    117.19;...
%    122.81;...
%    135.02;...
%    128.52];

% extract timecourse
fprintf('\nExtracting %g voxels from %g images\n',size(ROI.XYZvox,2),length(CurDes.V.vol))
for j = 1:length(CurDes.V.vol)
   if rem(j,100)==0
       fprintf('\nGot to image %g\n',j)
   end
   ROI.rawY(j,:) = spm_sample_vol(CurDes.V.vol(j), ROI.XYZvox(1,:), ROI.XYZvox(2,:), ROI.XYZvox(3,:),1);   
   % 1 at end of for trilinear interpolation - see
   % spm_sample_vol.
end

%average across voxels
ROI.rawY = mean(ROI.rawY,2);

% get baseline value
switch CurDes.base_flags.base_region
    case 1  % same ROI
        ROI.baseline = ROI.rawY;
        
    case 2  % grand mean of volume
        ROI.baseline = CurDes.V.G;  % global means
        
    case 3 % user-spec'ed mask
        % must extract values from mask
        baseV = spm_vol(CurDes.base_flags.base_mask);
        [baseY,base_XYZmm] = spm_read_vols(baseV);
        
        %select only voxels with value 1 in mask
        base_XYZmm = base_XYZmm(:,find(base_Y==max(max(max(baseY)))));
        
        %convert to vox
        vinv_base = inv(baseV.mat);  %matrix converting mm->voxels
        base_XYZvox = vinv_base(1:3,:)*[base_XYZmm; ones(1,size(base_XYZmm,2))];
        base_XYZvox = round(base_XYZvox);

        fprintf('\nExtracting %g baseline voxels from %g images\n',length(base_XYZvox),length(CurDes.V.vol))
        for j = 1:length(CurDes.V.vol)
            if rem(j,100)==0
                fprintf('\nGot to image %g\n',j)
            end
            ROI.baseline(j,1) = mean(spm_sample_vol(CurDes.V.vol(j), base_XYZvox(1,:), base_XYZvox(2,:), base_XYZvox(3,:),1));   
        end
        
    case 4 % user-spec'ed inverse mask
        % get info for inverse mask
        basenegV = spm_vol(CurrDes.base_flags.base_neg);
        [basenegY,baseneg_XYZmm] = spm_read_vols(basenegV);
        baseneg_XYZmm = baseneg_XYZmm(:,find(basenegY==max(max(max(basenegY)))));

        % get info for mask template
        basetempV = spm_vol(CurrDes.base_flags.base_template);
        [basetempY,basetemp_XYZmm] = spm_read_vols(basetempV);
        basetemp_XYZmm = basetemp_XYZmm(:,find(basetempY==max(max(max(basetempY)))));

        % select only voxels in template but not in inverse
        % mask.
        base_XYZmm = (setdiff(basetemplate_XYZmm', baseneg_XYZmm', 'rows'))';

        %convert to vox
        vinv_base = inv(basetempV.mat);
        base_XYZvox = vinv_base(1:3,:)*[base_XYZmm; ones(1, size(base_XYZmm, 2))];
        base_XYZvox = round(base_XYZvox);

        fprintf('\nExtracting %g baseline voxels from %g images\n',length(base_XYZvox),length(CurDes.V.vol))
        for j = 1:length(CurDes.V.vol)
            if rem(j,100)==0
                fprintf('\nGot to image %g\n',j)
            end
            ROI.baseline(j,1) = mean(spm_sample_vol(CurDes.V.vol(j), base_XYZvox(1,:), base_XYZvox(2,:), base_XYZvox(3,:),1));   
        end
        
end % switch over base_region

% that sets our baseline correctly for time-varying, but for
% steady, we must now average the baseline over time.

if ~CurDes.base_flags.base_timeflag
    % steady baseline
    if CurDes.base_flags.base_condflag
        % mean of condition - find which one and where it's
        % located in Cond vector - will be at several points
        % if we're allowing session-by-session.
        con_list = get(findobj('Tag', 'base_condflag_popup'), 'String');
        CurDes.base_flags.cond_name = con_list{CurDes.base_flags.base_cond_idx};
        [full_cond_names{1:length(CurDes.Cond)}] = deal(CurDes.Cond.name);
        cond_idx = find(strcmp(full_cond_names, CurDes.base_flags.cond_name));
        
        % accumulate stim vectors across sessions
        base_stim = zeros(length(CurDes.V.vol), 1);
        for c = 1:length(cond_idx)
            % we'll assume sessions don't overlap each
            % other, so each session's condition shouldn't
            % interfere with the others'.
            base_stim = base_stim + CurDes.Cond(cond_idx(c)).stim;
        end
        
        % filter basic baseline with aggregate stimulus
        % function from this condition.
        cond_baseline = ROI.baseline .* base_stim; 
        ROI.baseline = mean(cond_baseline);
    else
        % mean of timecourse
        ROI.baseline = mean(ROI.baseline);
    end
end


        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3.
% Convert into percent signal change and filter as
% appropriate (specified in interface structure).
% -----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert to % signal change.
% ---------------------------------

% If ROI.baseline is scalar, we can use it in array
% division, so this covers all our bases. (bases!  ha!)
ROI.rawY_perc = ROI.rawY ./ ROI.baseline;
ROI.rawY_perc = ROI.rawY_perc * 100;
ROI.rawY_perc = ROI.rawY_perc - 100;


% Filtering.
% -------------------

ROI.Y = ROI.rawY_perc;

% Detrending takes place within the deconvolution model, so
% we don't need to do it here (detrending will include a
% linear term in our model).

% Trimming first.
if CurDes.preproc_flags.trim
	 prev_idx = 0;
	 for s = 1:CurDes.nsess
		sess_idx = find(CurDes.V.iS==s);
		out_mask = (abs(zscore(ROI.Y(sess_idx)))) > CurDes.preproc_flags.trim_amt;
        out_idx = find(out_mask);
        in_idx = find(~out_mask);
        % don't want to include outliers in mean value we'll
        % replace them with, so we ignore them in mean.
		new_val = mean(ROI.Y(prev_idx+in_idx));
		ROI.Y(prev_idx+out_idx) = new_val;
		prev_idx = sess_idx(end);
	 end
 end
 
 
 % High-pass and low-pass next.
 % SPM2 and SPM99 have slightly different syntax for
 % spm_filter, so we have to take that into account.  Also,
 % SPM2 doesn't have low-pass filtering; we account for that
 % in roi_deconvolve_interface, but just in case...
 
 if CurDes.preproc_flags.lowpass | CurDes.preproc_flags.highpass
    for s = 1:CurDes.nsess
	    row{s} = find(CurDes.V.iS==s);
        K{s} = struct(	'HChoice',	CurDes.preproc_flags.high_choice,...
			'HParam',	CurDes.preproc_flags.high_param,...
			'LChoice',	CurDes.preproc_flags.low_choice,...
			'LParam',	(CurDes.preproc_flags.low_param)/sqrt(8*log(2)),...
			'row',		row{s},...
			'RT',		CurDes.TR);
    end
    filter_ver = which('spm_filter');
    if ~isempty(strfind(filter_ver, '99'));
        %i.e., we're in SPM99 right now
            K = spm_filter('set',K);    
    else  %i.e., SPM2
        % SPM2 filters aren't cell arrays - obviously.  I
        % mean, why be consistent between versions?  It's
        % not like anyone else will ever look at this code.
        % Oh, wait.
        for k = 1:length(K)
            newK(k) = K{k};
        end
        K = spm_filter(newK);
    end
end
if CurDes.preproc_flags.lowpass
    if ~isempty(strfind(filter_ver, '99'))
            %i.e., we're in SPM99
            ROI.Y = spm_filter('low',K,ROI.Y);
    else  %i.e., SPM2
            error('Low-pass filtering is no longer permitted in SPM2 - you shouldn''t have gotten here at all.');
    end
end
if CurDes.preproc_flags.highpass
    if ~isempty(strfind(filter_ver, '99'))
        %i.e., we're in SPM99
            ROI.Y = spm_filter('high',K,ROI.Y);
    else  %i.e., SPM2
            ROI.Y = spm_filter(K, ROI.Y);
    end
end

clear K;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 4.
% Deconvolve filtered timecourse into condition-by-condition
% impulse response functions (IRFs).
% -----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables
% ---------------------
n = length(CurDes.V.vol);
nsess = CurDes.nsess;
n_lags = CurDes.n_lags;
p = nsess; % one block parameter for each session
n_basic = nsess; % number of non-Cond params
if (CurDes.preproc_flags.detr)
    p = p + nsess; % one linear trend for each session
    n_basic = n_basic + nsess;
end

% n_lags parameters for each condition
p = p + (n_lags*length(CurDes.Cond));


% create X matrix of parameters - n_vols rows by n_params
% columns
%   - create block effects first - one for each session
%   - next linear effects if detr is on - just linearly
%   increasing numbers
%   - next loop over Cond structs
%       - n_lags columns for each
%   - then censor rows of X matrix from omitted

% Create X matrix
% ------------------------

X = zeros(n,p);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New Version Test here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Loop over Conds and add n_lags
%% time-shifted columns for each.
%for c = 1:length(CurDes.Cond)
%    for g = 1:n_lags
%        curr_idx = n_lags*(c-1) + g;
%        curr_stim = CurDes.Cond(c).stim;        
%        % add time-lag to top of stim fn
%        new_stim = [zeros(g-1,1); curr_stim];
%        % remove excess length
%        new_stim = new_stim(1:n);
%        % add to X matrix
%        X(:,n_basic + curr_idx) = new_stim;
%    end
%end

stim = horzcat(CurDes.Cond.stim);

for c = 1:size(stim,2)
    horz_offset = (n_basic + 1) + (n_lags * (c-1));
    X(:,horz_offset:(horz_offset + (n_lags - 1))) = repmat(stim(:,c),1,n_lags);
end

clear stim;

for j = (n_basic + 1):size(X,2)
    temp = X(:,j);
    offset = (rem((j - n_basic), n_lags) - 1);
    if offset < 0
        offset = n_lags - 1;
    end
    temp = [zeros(offset,1); temp];
    temp = temp(1:n);
    X(:,j) = temp;
end

% Censoring has been moved to roi_deconvolve_interface, into
% when stim functions are assembled.

%% Censor X matrix and V and Y vectors with it
%% -------------------------
%for s = 1:nsess
%    sess_idx = find(CurDes.V.iS==s);
%    % consistency check
%    if isempty(sess_idx)
%        error('Can''t find any rows for this session!');
%    end
%    % consistency check    
%    if (size(X,1)~=length(CurDes.V.vol))
%        error('X matrix and V vector are different sizes!');
%    end
%    first_sess_idx = sess_idx(1);
%    sess_offset = first_sess_idx - 1;
%    if ~isempty(CurDes.omitted)
%        censor = CurDes.omitted{s} + sess_offset;
%
%        CurDes.V.vol(censor) = [];
%        CurDes.V.iS(censor) = [];
%        CurDes.V.G(censor) = [];
%        ROI.Y(censor,:) = [];
%        
%        X(censor,:) = [];        
%    end
%    n = length(CurDes.V.vol);
%end

% add baseline effects finally
for s = 1:nsess
    sess_row = find(CurDes.V.iS==s);
    if CurDes.preproc_flags.detr
        sess_col = 1 + ((s - 1) * 2); % offset for linear effects
    else
        sess_col = s;
    end
    
    % first block effects
    X(sess_row,sess_col) = 1;
    
    % next linear effects if desired
    if CurDes.preproc_flags.detr
        X(sess_row,sess_col + 1) = (0:(length(sess_row) - 1))';
    end
end

disp('Created X matrix');

% Cache X matrix for full, base and reduced models
% ------------------------------------------------

% Full first
ROI.X = sparse(X);
clear X;

if debug
    disp('ROI.X');
    ROI.X
end


% Baseline model next
ROI.base_X = ROI.X(:,1:n_basic); % non-Cond effects
for c = 1:length(ROI.Cond)
    % add in any baseline stimulus conditions
    if ROI.Cond(c).basemodel_flag
        cond_idx = 1:n_lags + (n_basic + (n_lags*(c-1)));
        ROI.base_X = horzcat(ROI.base_X, ROI.X(:,cond_idx));
    end
end
if debug
    disp('ROI.base_x');
    ROI.base_X
end


% Reduced model for each Cond next
for c = 1:length(ROI.Cond)
    if ROI.Cond(c).basemodel_flag
        ROI.Cond(c).redcd_X = [];
    else
        % indices in X for this condition's effects
        cond_idx = (1:n_lags) + (n_basic + (n_lags*(c-1)));
        % remove these condition's indices in reduced model
        redcd_idx = setdiff(1:p, cond_idx);
        ROI.Cond(c).redcd_X = ROI.X(:, redcd_idx);
    end
end

if debug
    disp('ROI.Cond.redcd_x');
    for c = 1:length(ROI.Cond)
        ROI.Cond(c).redcd_X
    end
end

% Finally, cache degrees of freedom
ROI.df = length(ROI.Y) - size(ROI.X, 2); % num of timepoints - num of parameters


% create necessary inverted matrices
%   - get out (X'X)^-1 and ((X'X)^-1)X' for each of:
%       - full model
%       - baseline model alone
%       - each reduced model - one for each Cond
%   - get out C((X'X)^-1)C' for each linear contrast, as
%   well as constant matrix - which is as follows:
%       - (C((X'X)^-1)C')^-1 = temp1
%       - (((X'X)^-1)C')*temp1 = temp2
%       - temp2*c = temp3
%       - identity(size((X'X)^-1), 1) = temp4
%       - temp4 - temp3 = constant matrix desired.  Okay?
%       Okay.



% Cache inverted matrices for all the models we need
% ---------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We need to get rid of these lines - they're huge
% performance sucks, both in CPU time and memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First full
% inversion proceeds much faster with full matrices than
% with sparse matrices.
ROI.full_xtxinv = inv(full((ROI.X)' * ROI.X));
ROI.full_xtxinv = sparse(ROI.full_xtxinv);
ROI.full_xtxinvxt = ROI.full_xtxinv * (ROI.X)';
disp('Inverted X matrix!');

% Next baseline
ROI.base_xtxinv = inv(full((ROI.base_X)' * ROI.base_X));
ROI.base_xtxinv = sparse(ROI.base_xtxinv);
ROI.base_xtxinvxt = ROI.base_xtxinv * (ROI.base_X)';

% Last reduced models
for c = 1:length(ROI.Cond)
    if isempty(ROI.Cond(c).redcd_X)
        ROI.Cond(c).redcd_xtxinv = [];
        ROI.Cond(c).redcd_xtxinvxt = [];
    else
        ROI.Cond(c).redcd_xtxinv = inv(full((ROI.Cond(c).redcd_X)' * ROI.Cond(c).redcd_X));
        ROI.Cond(c).redcd_xtxinv = sparse(ROI.Cond(c).redcd_xtxinv);
        ROI.Cond(c).redcd_xtxinvxt = ROI.Cond(c).redcd_xtxinv * (ROI.Cond(c).redcd_X)';
    end
    if rem(c,10)==0
        fprintf('\nInverted condition %g matrix', c);
    end
end

% Cache inverted and constant matrices for contrasts
% -------------------------------------------------

for c = 1:length(ROI.Contrasts)
    ROI.Contrasts(c).cxtxinvct = (ROI.Contrasts(c).c * ROI.fullxtxinv) * (ROI.Contrasts(c).c)';
    tmp1 = inv(ROI.Contrasts(c).cxtxinvct);
    tmp2 = (ROI.fullxtxinv * (ROI.Contrasts(c).c)') * tmp1;
    tmp3 = tmp2 * ROI.Contrasts(c).c;
    tmp4 = eye(size(ROI.fullxtxinv, 1));
    ROI.Contrasts(c).constant = tmp4 - tmp3;
end


% do regression analysis
%   - calc betas for baseline with least-squares:
%   ((X'X)^-1)X'*ROI.Y = betas.
%   - calc sse for baseline: 
%       - X*b = y-hat
%       - y - y-hat = e
%       - e*e = sse
%   - calc betas for full model as above
%   - calc see for full as above, but also get out y-hat and
%   e as official fitted time series and residual
%   timeseries.
%   - calc t-coefficients + con ints for each ith beta:
%       - df = n(images) - p(parameters)
%       - mse = sse / df
%       - var = mse * ((X'X)^-1)[i, i]
%       - std = sqrt(var)
%       - t-coeff = beta / std - threshold with max and neg
%       max t-coeff
%       - conint = CON_COEFF * std
%   - determine significane for total model - F and r^2 as
%   below.
%   - determine significance of each Cond:
%       - cache irf and conints for irf
%       - calc betas for reduced as above
%       - calc sse for reduced as above
%       - calc F for reduced:
%           - p = # of params in full
%           - q = # of params in reduced
%           - msreg = (sse_red - sse_full) / (p - q)
%           - msef = sse_full / (n(images) - p)
%           - F = msreg / msef - threshold with max and 0
%       - calc r^2 for reduced:
%           - r^2 = (sse_red - sse_full) / sse_red
%           - threshold with 0 and 1



% Do least-square analysis for baseline and full
% ------------------------------------------------

% we use regress here instead of actual matrix
% multiplication for optimization with MATLAB.

% potential other code
%[ROI.betas_baseline, foo, ROI.base_err] = regress(ROI.Y, ROI.base_X);
%ROI.base_Y_est = ROI.base_X * ROI.betas_baseline;
%ROI.base_sse = dot(ROI.base_err,ROI.base_err);

% yet more potential other code
ROI.betas_baseline = ROI.base_X \ ROI.Y;
ROI.base_Y_est = ROI.base_X * ROI.betas_baseline;
ROI.base_err = ROI.Y - ROI.base_Y_est;
ROI.base_sse = dot(ROI.base_err , ROI.base_err);

% baseline first
%ROI.betas_baseline = ROI.base_xtxinvxt * ROI.Y;
%ROI.base_Y_est = ROI.base_X * ROI.betas_baseline;
%ROI.base_err = ROI.Y - ROI.base_Y_est;
%ROI.base_sse = dot(ROI.base_err , ROI.base_err);

% potential other code
%[ROI.betas_full, foo, ROI.full_err] = regress(ROI.Y, ROI.X);
%ROI.full_Y_est = ROI.X * ROI.betas_full;
%ROI.full_sse = dot(ROI.full_err,ROI.full_err);

% yet more potential other code
ROI.betas_full = ROI.X \ ROI.Y;
ROI.full_Y_est = ROI.X * ROI.betas_full;
ROI.full_err = ROI.Y - ROI.full_Y_est;
ROI.full_sse = dot(ROI.full_err, ROI.full_err);

% Next full model
%ROI.betas_full = ROI.full_xtxinvxt * ROI.Y;
%ROI.full_Y_est = ROI.X * ROI.betas_full;
%ROI.full_err = ROI.Y - ROI.full_Y_est;
%ROI.full_sse = dot(ROI.full_err, ROI.full_err);

% Calculate t-stats and con ints for each beta
for b = 1:length(ROI.betas_full)
    curr_df = n - p;
    curr_mse = ROI.full_sse / curr_df;
    % is there a way we could do this without the inverse?
    curr_var = curr_mse * ROI.full_xtxinv(b,b);
    curr_std = sqrt(curr_var);
    if curr_std < tol
        if (ROI.betas_full(b) > 0) ROI.tstat(b) = Inf;
        else ROI.tstat(b) = -Inf; end
    else
        ROI.tstat(b) = ROI.betas_full(b) / curr_std;
    end
    ROI.conint(b) = CON_COEFF * curr_std;
    ROI.stderr(b) = curr_std;
end

% Calculate significance of total model - F and r^2
% --------------------------------------------------

% F first
p_reduced = size(ROI.base_X,2);
if p_reduced > p
    error('Baseline model larger than full model - check code.')
end
ms_regression = (ROI.base_sse - ROI.full_sse) / (p - p_reduced);
if (ms_regression < 0)   ms_regression = 0; end

ms_error_full = ROI.full_sse / (n - p);
if (ms_error_full < 0)  ms_error_full = 0; end

if ms_error_full < tol  
    ROI.F_total = Inf;
else  
    ROI.F_total = ms_regression / ms_error_full;
end
if (ROI.F_total < 0)    ROI.F_total = 0;    end

% Next r^2
if (ROI.base_sse < tol)  
    ROI.r_sqr_total = 0;
else    
    ROI.r_sqr_total = (ROI.base_sse - ROI.full_sse) / ROI.base_sse; 
end

if (ROI.r_sqr_total > 1)  ROI.r_sqr_total = 1;
elseif (ROI.r_sqr_total < 0)   ROI.r_sqr_total = 0;  end

% Calculate significance of each condition - F and r^2
% ----------------------------------------------------

for c = 1:length(ROI.Cond)
    
    % cache impulse response function and con ints first
    % Hey, look, these six lines are the point of this whole
    % program!  Awesome!
    cond_idx = (1:n_lags) + (n_basic + (n_lags * (c-1)));
    ROI.Cond(c).irf = ROI.betas_full(cond_idx);
    ROI.Cond(c).conint = ROI.conint(cond_idx);
    ROI.Cond(c).stderr = ROI.stderr(cond_idx);
    ROI.Cond(c).tstat = ROI.tstat(cond_idx);
    ROI.Cond(c).df = ROI.df;
    
    if ROI.Cond(c).basemodel_flag
        % condition is in baseline model, so we don't care
        % about its significance individually
        ROI.Cond(c).betas_redcd = [];
        ROI.Cond(c).redcd_Y_est = [];
        ROI.Cond(c).redcd_err = [];
        ROI.Cond(c).redcd_sse = [];
        ROI.F_redcd = [];
        ROI.r_sqr_redcd = [];
        continue;
    end
    
    % calc reduced betas next

    % potential other code
    [ROI.Cond(c).betas_redcd, foo, ROI.Cond(c).redcd_err] = regress(ROI.Y, ROI.Cond(c).redcd_X);
    ROI.Cond(c).redcd_Y_est = ROI.Cond(c).redcd_X * ROI.Cond(c).betas_redcd;
    ROI.Cond(c).redcd_sse = dot(ROI.Cond(c).redcd_err, ROI.Cond(c).redcd_err);
    
    % yet more potential other code
    ROI.Cond(c).betas_redcd = ROI.Cond(c).redcd_X \ ROI.Y;
    ROI.Cond(c).redcd_Y_est = ROI.Cond(c).redcd_X * ROI.Cond(c).betas_redcd;
    ROI.Cond(c).redcd_err = ROI.Y - ROI.Cond(c).redcd_Y_est;
    ROI.Cond(c).redcd_sse = dot(ROI.Cond(c).redcd_err , ROI.Cond(c).redcd_err);


    %ROI.Cond(c).betas_redcd = ROI.Cond(c).redcd_xtxinvxt * ROI.Y;
    %ROI.Cond(c).redcd_Y_est = ROI.Cond(c).redcd_X * ROI.Cond(c).betas_redcd;
    %ROI.Cond(c).redcd_err = ROI.Y - ROI.Cond(c).redcd_Y_est;
    %ROI.Cond(c).redcd_sse = dot(ROI.Cond(c).redcd_err , ROI.Cond(c).redcd_err);

    % now calc significance for condition
    p_reduced = size(ROI.Cond(c).redcd_X,2);
    if p_reduced > p
        error('Condition reduced model larger than full model - check code.')
    end
    ms_regression = (ROI.Cond(c).redcd_sse - ROI.full_sse) / (p - p_reduced);
    if (ms_regression < 0)   ms_regression = 0; end

    ms_error_full = ROI.full_sse / (n - p);
    if (ms_error_full < 0)  ms_error_full = 0; end

    if ms_error_full < tol  
        ROI.Cond(c).F_redcd = Inf;
    else  
        ROI.Cond(c).F_redcd = ms_regression / ms_error_full;
    end
    if (ROI.Cond(c).F_redcd < 0)    ROI.Cond(c).F_redcd = 0;    end

    % Finally, r^2
    if (ROI.Cond(c).redcd_sse < tol)  
        ROI.Cond(c).r_sqr_redcd = 0;
    else    
        ROI.Cond(c).r_sqr_redcd = (ROI.Cond(c).redcd_sse - ROI.full_sse) / ROI.Cond(c).redcd_sse; 
    end

    if (ROI.Cond(c).r_sqr_redcd > 1)  ROI.Cond(c).r_sqr_redcd = 1;
    elseif (ROI.Cond(c).r_sqr_redcd < 0)   ROI.Cond(c).r_sqr_redcd = 0;  end
    
    
end % condition loop

% do GLT analysis
%   for each desired GLT:
%   - calc linear combination - C * betas
%   - calc t-coeff for each linear comb - as above.
%   - calc r-coeff for reduced model:
%       - multiply constant matrix a above: a*betas =
%       r-coeff.
%   - calc sse_reduced on r-coeff.
%   - calc F for this reduced GLT - using sse_reduced from
%   r-coeffs - as above
%   - calc r^2 for this reduced GLT - using sse_reduced from
%   r-coeffs - as above 

% Calculate significance of each contrast - t, F and r^2
% ------------------------------------------------------

for c = 1:length(ROI.Contrasts)
    
    % potential other code - use above
    %[ROI.Contrasts(c).betas_redcd, foo, ROI.Contrasts(c).redcd_err] = regress(ROI.Y, ROI.Contrasts(c).redcd_X);
    %ROI.Contrasts(c).redcd_Y_est = ROI.Contrasts(c).redcd_X * ROI.Contrasts(c).betas_redcd;
    %ROI.Contrasts(c).redcd_sse = dot(ROI.Contrasts(c).redcd_err, ROI.Contrasts(c).redcd_err);
    

    % calc reduced betas first
    ROI.Contrasts(c).lin_comb = ROI.Contrasts(c).c * ROI.betas_full;
    % lin_comb has as many elements as c has rows
    ROI.Contrast(c).betas_con = ROI.Contrast(c).constant * ROI.Y;
    ROI.Contrast(c).con_Y_est = ROI.X * ROI.betas_con;
    ROI.Contrast(c).con_err = ROI.Y - ROI.Contrast(c).con_Y_est;
    ROI.Contrast(c).con_sse = dot(ROI.Contrast(c).con_err , ROI.Contrast(c).con_err);
    ROI.Contrast(c).df = ROI.df;

    % calc t-stats for each constraint / lin_comb element
    for t = 1:length(ROI.Contrast(c).lin_comb)
        curr_df = n - p;
        curr_mse = ROI.full_sse / curr_df;
        curr_var = curr_mse * ROI.Contrast(c).cxtxinvct(t,t);
        curr_std = sqrt(curr_var);
        if curr_std < tol
            if (ROI.Contrast(c).lin_comb(t) > 0) ROI.Contrast(c).lin_tstat(t) = Inf;
            else ROI.Contrast(c).lin_tstat(t) = -Inf; end
        else
            ROI.Contrast(c).lin_tstat(t) = ROI.Contrast(c).lin_comb(t) / curr_std;
        end
    end
    
    % now calc significance for contrast test
    p_reduced = p - length(ROI.Contrast(c).lin_comb);
    if p_reduced < p
        error('Reduced model larger than full model - check code.')
    end
    ms_regression = (ROI.Contrast(c).con_sse - ROI.full_sse) / (p - p_reduced);
    if (ms_regression < 0)   ms_regression = 0; end

    ms_error_full = ROI.full_sse / (n - p);
    if (ms_error_full < 0)  ms_error_full = 0; end

    if ms_error_full < tol  
        ROI.Contrast(c).F_con = Inf;
    else  
        ROI.Contrast(c).F_con = ms_regression / ms_error_full;
    end
    if (ROI.Contrast(c).F_con < 0)    ROI.Contrast(c).F_con = 0;    end

    % Next r^2
    if (ROI.Contrast(c).con_sse < tol)  
        ROI.Contrast(c).r_sqr_con = 0;
    else    
        ROI.Contrast(c).r_sqr_con = (ROI.Contrast(c).con_sse - ROI.full_sse) / ROI.Contrast(c).con_sse; 
    end

    if (ROI.Contrast(c).r_sqr_con > 1)  ROI.Contrast(c).r_sqr_con = 1;
    elseif (ROI.Contrast(c).r_sqr_con < 0)   ROI.Contrast(c).r_sqr_con = 0;  end
    
end % contrast loop


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 5.
% Create output files for this analysis.
% -----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Output will be :
%   - one .txt file containing whole filtered percent signal timecourse.
%   - one .txt file containing information for each
%   condition, namely the following:
%      - separate impulse response functions for each condition.
%       - std deviations for each timepoint of the IRF.
%       - conint for each timepoint of the IRF
%       - t-stats and p-stats for each timepoint different from
%       the baseline model.
%       - F-stats and r^2 for each condition against baseline
%       model.
%       - F-stat and r^2 for each GLT (as well as t-stat for
%       each constraint in the GLT).
%   - one .mat file containing the ROI analysis structure
%   (or possibly some summary structure, if it ends up being
%   huge).
% 
% Note that each Cond structure is a condition in one
% session, so in a multi-session experiment, we must
% accumulate all the info over all sessions.
%
% If requested, we can produce a separate
% condition-by-condition file which breaks out the
% conditions by session.

cdir = pwd;
cd(CurDes.outdir);
% ROI.name is complete subject x ROI plus prepend string
filename = [ROI.name '_percent_signal'];
filename_cond = [filename '_condition'];
filename_stats = [filename '_stats'];

% Save whole timecourse
% ---------------------
fid = fopen([filename '.txt'],'w');
if fid == -1
    error(['Unable to create new file - please check permissions in current directory.']); 
end
fprintf(fid,'%4.4f\n', ROI.Y);
fclose(fid);
   
% Accumulate conditions over sessions
% ---------------------------------------
clear cond_idx;
[cond_idx{1:length(ROI.Cond)}] = deal(ROI.Cond.sess);

% this is a dummy struct, to allow the code below
avg_Cond = struct('name', [], 'irf', [], 'conint', [], 'stderr', [], 'tstat', [], 'df', [], 'basemodel_flag', 0);
for c = 1:length(ROI.Cond)
    avg_idx = 0;
    for a = 1:length(avg_Cond)
        if strcmp(avg_Cond(a).name, ROI.Cond(c).name)
            avg_idx = a;
        end
    end
    if avg_idx==0
        % no avg_Cond for this condition exists yet, so
        % let's make one
        avg_Cond(end + 1) = struct('name', ROI.Cond(c).name, ...
            'irf', [], ...
            'conint', [], ...
            'stderr', [], ...
            'tstat', [], ...
            'df', [], ...
            'basemodel_flag', 0);
        avg_idx = length(avg_Cond);
    end
    % tack on this condition's stuff to its group
    % structure
    avg_Cond(avg_idx).irf = horzcat(avg_Cond(avg_idx).irf, ROI.Cond(c).irf);
    avg_Cond(avg_idx).conint = horzcat(avg_Cond(avg_idx).conint, ROI.Cond(c).conint');
    avg_Cond(avg_idx).stderr = horzcat(avg_Cond(avg_idx).stderr, ROI.Cond(c).stderr');
    avg_Cond(avg_idx).tstat = horzcat(avg_Cond(avg_idx).tstat, ROI.Cond(c).tstat');
end
% strip dummy struct
avg_Cond(1) = [];

% now avg_Cond is a list of all the named conditions, with a
% column in each of its stats structures corresponding to
% each session that a Cond with its name occurred in.
%   So all we need do is accumulate across columns to
% accumulate across sessions...
%   We could do some kind of trickier combination of these
% than just a simple average - a repeated-measures statistic
% or something? - but for now the simple average will do.

for a = 1:length(avg_Cond)
    avg_Cond(a).irf = mean(avg_Cond(a).irf, 2);
    avg_Cond(a).conint = mean(avg_Cond(a).conint, 2);
    avg_Cond(a).stderr = mean(avg_Cond(a).stderr, 2);
    avg_Cond(a).tstat = mean(avg_Cond(a).tstat, 2);
    avg_Cond(a).df = ROI.df;
end
ROI.avg_Cond = avg_Cond;

% Save condition-by-condition information
% ---------------------------------------
fid = fopen([filename_cond '.txt'],'w');
if fid == -1
    error(['Unable to create new file - please check permissions in current directory.']); 
end

% Baseline first - same for all conditions
fprintf(fid, 'Baseline:  ');
if CurDes.base_flags.base_timeflag
    fprintf(fid, 'time-varying,  ');
else 
    fprintf(fid, 'steady,  ');
end
if CurDes.base_flags.base_condflag
    fprintf(fid, 'mean of condition %s,  ', CurDes.base_flags.cond_name);
else 
    fprintf(fid, 'mean of timecourse,  ');
end
if CurDes.base_flags.base_region==1
    fprintf(fid, 'same ROI.\n\n');
elseif CurDes.base_flags.base_region==2
    fprintf(fid, 'grand mean of volume.\n\n');
elseif CurDes.base_flags.base_region==3
    fprintf(fid, 'mean of user-specified mask.\n');
    fprintf(fid, 'Baseline mask: %s\n\n', CurDes.base_flags.base_mask);
else
    fprintf(fid, 'mean of user-specified inverse mask.\n');
    fprintf(fid, 'Baseline template (included vox): %s\n', CurDes.base_flags.base_template);
    fprintf(fid, 'Baseline negative mask (excluded mask): %s\n\n', CurDes.base_flags.base_neg);
end

% next each condition
for c = 1:length(ROI.avg_Cond)
    fprintf(fid,'Condition:\t%s', ROI.avg_Cond(c).name);
    fprintf(fid,'\nIRF (%% sig. change):\t');
    fprintf(fid,'%4.4f\t', ROI.avg_Cond(c).irf);
    fprintf(fid,'\n%g%% confidence int. 1/2 width:\t', CON_INT);
    fprintf(fid,'%4.4f\t', ROI.avg_Cond(c).conint);
    fprintf(fid,'\nstd. error:\t');
    fprintf(fid,'%4.4f\t', ROI.avg_Cond(c).stderr);
    fprintf(fid,'\nt-stat (diff. from baseline model):\t');
    fprintf(fid,'%4.4f\t', ROI.avg_Cond(c).tstat);
    fprintf(fid,'\n\n');
    % save reduced model fit and r^2?
end
fclose(fid);



% Save session-by-session condition info (if requested)
% ------------------------------------------------------
if CurDes.sess_out

fid = fopen([filename_cond '_sess.txt'],'w');
if fid == -1
    error(['Unable to create new file - please check permissions in current directory.']); 
end

% Baseline first - same for all conditions
fprintf(fid, 'Baseline:  ');
if CurDes.base_flags.base_timeflag
    fprintf(fid, 'time-varying,  ');
else 
    fprintf(fid, 'steady,  ');
end
if CurDes.base_flags.base_condflag
    fprintf(fid, 'mean of condition %s,  ', CurDes.base_flags.cond_name);
else 
    fprintf(fid, 'mean of timecourse,  ');
end
if CurDes.base_flags.base_region==1
    fprintf(fid, 'same ROI.\n\n');
elseif CurDes.base_flags.base_region==2
    fprintf(fid, 'grand mean of volume.\n\n');
elseif CurDes.base_flags.base_region==3
    fprintf(fid, 'mean of user-specified mask.\n');
    fprintf(fid, 'Baseline mask: %s\n\n', CurDes.base_flags.base_mask);
else
    fprintf(fid, 'mean of user-specified inverse mask.\n');
    fprintf(fid, 'Baseline template (included vox): %s\n', CurDes.base_flags.base_template);
    fprintf(fid, 'Baseline negative mask (excluded mask): %s\n\n', CurDes.base_flags.base_neg);
end

% next each condition
for c = 1:length(ROI.Cond)
    fprintf(fid,'Condition:\t%s,session %g', ROI.Cond(c).name, ROI.Cond(c).sess);
    fprintf(fid,'\nIRF (%% sig. change):\t');
    fprintf(fid,'%4.4f\t', ROI.Cond(c).irf);
    fprintf(fid,'\n%g%% confidence int. 1/2 width:\t', CON_INT);
    fprintf(fid,'%4.4f\t', ROI.Cond(c).conint);
    fprintf(fid,'\nstd. error:\t');
    fprintf(fid,'%4.4f\t', ROI.Cond(c).stderr);
    fprintf(fid,'\nt-stat (diff. from baseline model):\t');
    fprintf(fid,'%4.4f\t', ROI.Cond(c).tstat);
    fprintf(fid,'\n\n');
    % save reduced model fit and r^2?
end
fclose(fid);

end % if CurDes.sess_out

% Save full-model and linear contrast information
% ------------------------------------------------
fid = fopen([filename_stats '.txt'],'w');
if fid == -1
    error(['Unable to create new file - please check permissions in current directory.']); 
end

% Full-model stuff first
fprintf(fid, 'Full model information\n');
fprintf(fid, '-----------------------\n');
fprintf(fid, 'Full model fit (F): %4.4f\n', ROI.F_total);
fprintf(fid, 'Full model fit (r-squared): %4.4f\n\n', ROI.r_sqr_total);

% Next contrast info
for c = 1:length(ROI.Contrasts)
    fprintf(fid,'Contrast:\t%g', c);
    temp_c = full(ROI.Contrasts(c).c);
    for n = 1:size(ROI.Contrasts(c).c,1)
        fprintf(fid,'\nContrast matrix:\t%g\t', temp_c(n,:));
    end
    clear temp_c;
    fprintf(fid,'\nContrast t-stats:\t');
    fprintf(fid,'%4.4f\t', ROI.Contrasts(c).lin_tstat);
    fprintf(fid,'\nContrast fit (F):\t%4.4f', ROI.Contrasts(c).F_con);
    fprintf(fid,'\nContrast fit (r-squared):\t%4.4f', ROI.Contrasts(c).r_sqr_con);
    fprintf(fid,'\n\n');
end


% Save X matrix, for reference, and/or use in AFNI
% comparison
temp_x = full(ROI.X);
fprintf(fid, 'X Matrix:\n');
for x = 1:size(ROI.X, 1)
    fprintf(fid, '%g\t', temp_x(x,:));
    fprintf(fid, '\n');
end
fprintf(fid, '\n\n');
fclose(fid);
clear temp_x;

% Save output structures in .mat file and give feedback
% ------------------------------------------------------

save([filename '.mat'], 'ROI', '-mat');
fprintf('\nThe percent signal change vector was saved as \n%s.mat \nand %s.txt.\n\n',...
    fullfile(pwd,filename), fullfile(pwd,filename));


% accumulate ROI struct, for averaging over group
% ------------------------------------------------

% make this an array, with a row for each ROI, and a column
% for each subject

vec_idx = 0;
for r = 1:size(ROI_vector, 1)
    if strcmp(ROI_vector{r}(1).roi_name, ROI.roi_name)
        % row exists for this ROI - tack on to end
        vec_idx = r;
        ROI_vector{vec_idx}(end + 1) = ROI;
    end
end
if vec_idx==0
    % no row yet for this ROI - make one
    ROI_vector = vertcat(ROI_vector, cell(1));
    vec_idx = length(ROI_vector);
    ROI_vector{vec_idx} = ROI;
end


end % Design struct (subject x ROI) loop


%%%%%%%%%%%%%%%%%%
% Step 6 - Group averaging - missing here until next week,
% when I'll attempt to finish it.  It's mostly written, but
% badly so far.
%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 7.
% Launch roi_deconvolve viewer
% -----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check first to see if viewer is already launched - don't
% re-launch if one is already up.

% this bit disabled till Group is back in.
%if isempty(findobj('Tag', 'Viewer'))
%    vec = ROI_vector(1);  % pick first ROI arbitrarily
%    if length(vec) < 2
%        % single-subject output only
%        fname = [vec{1}.name '_percent_signal.mat'];
%    else
%        % group output available
%        fname = [vec{1}.roi_name '_percent_signal_group.mat'];
%    end
%    
%    roi_deconvolve_viewer('init', fname);
%end

if isempty(findobj('Tag', 'Viewer'))
    roi_deconvolve_viewer('init', [ROI.name '_percent_signal.mat']);
end


% return to cached present directory
cd(cdir);

return;