%RRI_ANALYSIS_BOOT Apply Behavioral or Task PLS test and bootstrap test
%	on RRI scan.
%
%   Usage: [brainlv,s,behavlv,brainscores,behavscores,lvcorrs,boot_result]= ...
%                rri_analysis_boot(ismean, ishelmert, iscontrast, isbehav, ...
%                newdata_lst,num_cond_lst,num_subj_lst, ...
%                behavdata_lst, helmertdata_lst, contrastdata_lst, ...
%                num_boot, Clim);
%
%   See also PLS_BOOT_TEST, PLS_DEVIATION_BOOT_TEST, BEHAVPLS_BOOT, TASKPLS_BOOT
%

%   I (ismean) - 1 if using grand mean deviation;
%   I (ishelmert) - 1 if using helmert matrix;
%   I (iscontrast) - 1 if using contrast data;
%   I (isbehav) - 1 if using behavior data;
%   I (newdata_lst) - A group list of datamat files;
%   I (num_cond_lst) - A group list of condition numbers;
%   I (num_subj_lst) - A group list of subject numbers;
%   I (behavdata_lst) - A group list of behav data with selected columns;
%   I (helmertdata_lst) - A group list of helmert matrix with selected columns;
%   I (contrastdata_lst) - A group list of contrast data with selected columns;
%                       contrastdata_lst has been orthchecked in
%                       erp_get_common function
%   I (num_boot) - Number of Permutation;
%   I (Clim) - upper limit of confidence interval estimated.
%
%   O (brainlv) - Left singular value vector. It is LV for the brain dimension.
%   O (s) - Singular values vector.
%   O (behavlv) - Right singular vector. It is LV for the behavior OR contrast
%		dimension for each scan.
%   O (brainscores) - The brain score.
%   O (behavscores) - If use GRAND_MEAN method, behavscores = behavlv, then
%		expand each condition for all the subjects. If not, it is
%		calculated for each set of condition, then stack together.
%   O (lvcorrs) - Correlates brain scores with behavior data for multiple scans,
%		return [] if for Task PLS or if use GRAND_MEAN method.
%   O (boot_result) - A Structure array containing the bootstrap result data.
%
%   Created on 03-OCT-2002 by Jimmy Shen for PLS test
%   Modified on 27-OCT-2002 by Jimmy Shen to add bootstrap test
%   Modifyed on Jan 13,03 for ERP analysis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs,boot_result,datamatcorrs_lst] = ...
                rri_analysis_boot(ismean, ishelmert, iscontrast, isbehav, ...
                newdata_lst,num_cond_lst,num_subj_lst, ...
                behavdata_lst, helmertdata_lst, contrastdata_lst, ...
                num_boot, Clim);

    % Init
    %
    brainlv = [];
    s = [];
    behavlv = [];
    brainscores = [];
    behavscores = [];
    lvcorrs = [];
    boot_result = [];

    datamatcorrs_lst = {};

    stacked_behavdata = [];
    stacked_helmertdata = [];
    stacked_contrastdata = [];
    stacked_datamatcorrs = [];
    stacked_datamat = [];
    stacked_mask = [];

    num_groups = length(newdata_lst);

    % keeps track of number of times a new bootstrap had to be generated
    %
    countnewboot=0;
    badbeh = [];

    progress_hdl = rri_progress_ui('initialize');

    msg = 'Working on PLS ...';
    rri_progress_ui(progress_hdl, '', msg);

    boot_progress = rri_progress_ui('initialize');

    if ismean
       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1, boot_progress);
    else
       [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 0, boot_progress);
    end

    if isempty(reorder)
       return;
    end;

    % loop accross the groups, and
    % calculate datamatcorrs for each group
    %
    for i = 1:num_groups

        k = num_cond_lst(i);
        n = num_subj_lst(i);

        datamat = newdata_lst{i};

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+1/(10*num_groups));

        % compute correlation or covariance
        %
        if ismean
            datamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat);
        elseif ishelmert
            datamatcorrs = rri_xcovy(helmertdata_lst{i}, datamat);
        elseif iscontrast
            datamatcorrs = rri_xcovy(contrastdata_lst{i}, datamat);
        elseif isbehav
            datamatcorrs = rri_corr_maps(behavdata_lst{i}, datamat, n, k);
            datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}];
        end

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+3/(10*num_groups));

	% if more than one group, stack data together
	%
        if ishelmert
            stacked_helmertdata = [stacked_helmertdata; helmertdata_lst{i}];
        elseif iscontrast
            stacked_contrastdata=[stacked_contrastdata;contrastdata_lst{i}];
        elseif isbehav
            stacked_behavdata = [stacked_behavdata; behavdata_lst{i}];
        end
        stacked_datamat = [stacked_datamat; datamat];
        stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];

        rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+5/(10*num_groups));

    end		% for

    % actually, all the groups must have the same condition number
    %
    num_cond = num_cond_lst(1);

    % Singular Value Decomposition
    %
    [r c] = size(stacked_datamatcorrs);
    if r <= c
        % transpose datamatcorrs to ensure SVD operation will be
        % on smallest of RxC dimension
        %
        [brainlv,s,behavlv] = svd(stacked_datamatcorrs',0);
    else
        [brainlv,s,behavlv] = svd(stacked_datamatcorrs,0);
    end

    original_sal = brainlv * s;
    orig_behavlv = behavlv * s;

    s = diag(s);

    rri_progress_ui(progress_hdl,'',9/10);

    % calculate behav scores
    %
    if ismean

        brainscores = stacked_datamat * brainlv;

        % Here, behavlv is actually designlv
        % according to taskpls.m: fscores=design=testvec(designlv)
        % so, designscores = designlv

        num_col = size(behavlv, 2);

        % expand the num_subj for each row (cond)
        % did the samething as testvec
        %
        for i = 1:num_groups

            k = num_cond_lst(i);
            n = num_subj_lst(i);

            tmp = reshape(behavlv((i-1)*k+1:(i-1)*k+k,:), ...
                [1, num_col*k]);
            tmp = repmat(tmp, [n, 1]);          % expand to num_subj
            tmp = reshape(tmp, [n*k, num_col]);

            behavscores = [behavscores; tmp];   % stack by groups

        end

    elseif ishelmert | iscontrast
        if ishelmert
            behav = stacked_helmertdata;
        elseif iscontrast
            behav = stacked_contrastdata;
        end

        [brainscores, behavscores] = ...
                rri_get_designscores(stacked_datamat, behav, ...
                brainlv, behavlv, num_cond_lst(1), num_subj_lst);
    elseif isbehav
        behav = stacked_behavdata;

        [brainscores, behavscores, lvcorrs] = ...
		rri_get_behavscores(stacked_datamat, behav, ...
		brainlv, behavlv, num_cond_lst(1), num_subj_lst);

        orig_corr = lvcorrs;
        [r1 c1] = size(orig_corr);
        distrib = zeros(r1, c1, num_boot+1);
        distrib(1:r1, 1:c1, 1) = orig_corr;
    end         % if

    rri_progress_ui(progress_hdl,'',1);

    %  Begin bootstrap loop
    %
    %  reorder = rri_mkboot_order(num_cond_lst(1), num_subj_lst, num_boot);

    %  move the following code up, so we don't need to run obserbation
    %  test if bootstrap can not run
    %
    % [reorder, new_num_boot] = ...
	% rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1);

    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 isbehav
            distrib = zeros(r1, c1, num_boot+1);
            distrib(1:r1, 1:c1, 1) = orig_corr;
        end
    end

    max_subj_per_group = 8;
    if (sum(num_subj_lst <= max_subj_per_group) == num_groups)
        is_boot_samples = 1;
    else
        is_boot_samples = 0;
    end

    if isempty(reorder)
       return;
    end

    if ~isbehav
        sal_sq = zeros(size(brainlv));
        sal_sum = zeros(size(brainlv));
        dsal_sq = zeros(size(behavlv));
        dsal_sum = zeros(size(behavlv));
    else
        sal_sq = original_sal.^2;
        sal_sum = original_sal;
        dsal_sq = orig_behavlv.^2;
        dsal_sum = orig_behavlv;
    end

    for p=1:num_boot

        msg = ['Working on Bootstrap:  ',num2str(p),' out of ',num2str(num_boot)];
        rri_progress_ui(progress_hdl, '', msg);
        rri_progress_ui(progress_hdl,'',p/num_boot);

        data_p = stacked_datamat(reorder(:,p),:);
        if ishelmert
            behav_p = stacked_helmertdata(reorder(:,p),:);
        elseif iscontrast
            behav_p = stacked_contrastdata(reorder(:,p),:);
        elseif isbehav
            % behav_p = stacked_behavdata(reorder(:,p),:);
        end
        stacked_data = [];

        for g=1:num_groups

            k = num_cond_lst(g);
            n = num_subj_lst(g);
            span = sum(num_subj_lst(1:g-1)) * num_cond;	% group length

            if ismean
                if num_groups == 1
                   meanmat = rri_task_mean(data_p, n);
                   data = meanmat - (ones(k,1)*mean(meanmat));
                else
%                   data = rri_task_mean(data_p(1+span:n*k+span,:), n) ...
%			- ones(k,1)*mean(data_p(1+span:n*k+span,:));
                   meanmat = rri_task_mean(data_p(1+span:n*k+span,:), n);
                   data = meanmat - (ones(k,1)*mean(meanmat));
                end
            elseif ishelmert
                if num_groups == 1
                   data = rri_xcovy(behav_p, data_p);
                else
                   data = rri_xcovy(behav_p(1+span:n*k+span,:), ...
				data_p(1+span:n*k+span,:));
                end
            elseif iscontrast
                if num_groups == 1
                   data = rri_xcovy(behav_p, data_p);
                else
                   data = rri_xcovy(behav_p(1+span:n*k+span,:), ...
				data_p(1+span:n*k+span,:));
                end
            elseif isbehav

                if ~is_boot_samples

                   % the code below is mainly trying to find a proper
                   % reorder matrix

                   % init badbehav cell array to 0
                   % which is used to record the bad behav data caused by
                   % bad re-order. This var. is for disp only.
                   %
                   badbehav = zeros(num_cond, size(behavdata_lst{g},2));

                   % Check for upcoming NaN and re-sample if necessary.
                   % this only happened on behavior analysis, because the
                   % 'xcor' inside of 'corr_maps' contains a 'stdev', which
                   % is a divident. If it is 0, it will cause divided by 0
                   % problem.
                   % since this happend very rarely, so the speed will not
                   % be affected that much.
                   %
                   % For behavpls_boot, also need to account for multiple
                   % scans and behavs
                   %
                   for c=1:k	% traverse all conditions in this group
                       stdmat(c,:) = std(stacked_behavdata(reorder((1+ ...
                          (n*(c-1))+span):(n*c+span), p), :));
                   end		% scanloop

                   % now, check to see if any are zero
                   %
                   while sum(stdmat(:)==0)>0
                       countnewboot = countnewboot + 1;

                       % keep track of scan & behav that force a resample
                       %
                       badbehav(find(stdmat(:)==0)) = ...
                          badbehav(find(stdmat(:)==0)) + 1;

                       % num_boot is just something to be picked to prevent
                       % infinite loop
                       %
                       if countnewboot > num_boot
                           msgbox(['countnewboot exceeds num_boot: ', ...
                              num2str(countnewboot)],'modal');
                           breakon=1;
                           break;
                       end

                       % reorder(:,p) = rri_mkboot_order(num_cond_lst(1), ...
                       %    num_subj_lst, 1);

                       boot_progress = rri_progress_ui('initialize');
                       reorder(:,p) = rri_boot_order(num_subj_lst, num_cond_lst(1), 1, 0, boot_progress);

                       for c=1:k	% recalc stdmat
                           stdmat(c,:) = std(stacked_behavdata(reorder((1+ ...
                              (n*(c-1))+span):(n*c+span), p), :));
                       end	% scanloop

                   end		% while

                   badbeh{g} = badbehav;	% save instead of disp

                end	% if ~is_boot_samples

		% now, we can use this proper reorder matrix
		% to generate behav_p & data_p, and then
		% to calculate datamatcoors
		%
                behav_p = stacked_behavdata(reorder(:,p),:);
                data_p = stacked_datamat(reorder(:,p),:);

                if num_groups == 1
                   data = rri_corr_maps(behav_p, data_p, n, k);
                else
                   data = rri_corr_maps(behav_p(1+span:n*k+span,:), ...
				data_p(1+span:n*k+span,:), n, k);
                end
            end		% ismean,isbehav, ...

            stacked_data = [stacked_data; data];

        end	% for num_groups

        % Singular Value Decomposition
        %
        [r c] = size(stacked_data);
        if r <= c
            % transpose datamatcorrs to ensure SVD operation will be
            % on smallest of RxC dimension
            %
            [pbrainlv, sboot, pbehavlv] = svd(stacked_data',0);
        else
            [pbrainlv, sboot, pbehavlv] = svd(stacked_data,0);
        end

        %  rotate pbehavlv to align with the original behavlv
        %
        rotatemat = rri_bootprocrust(behavlv,pbehavlv);

        %  rescale the bootstrap vectors
        %
        pbrainlv = pbrainlv * sboot * rotatemat;
        pbehavlv = pbehavlv * sboot * rotatemat;

        if isbehav
            [brainsctmp, behavsctmp, bcorr] = ...
		rri_get_behavscores(data_p, behav_p, rri_normalize(pbrainlv), ...
		rri_normalize(pbehavlv), num_cond_lst(1), num_subj_lst);

            distrib(1:r1, 1:c1, p+1) = bcorr;
        end

        sal_sum = sal_sum + pbrainlv;
        sal_sq = sal_sq + pbrainlv.^2;

        dsal_sum = dsal_sum + pbehavlv;
        dsal_sq = dsal_sq + pbehavlv.^2;

    end		% for num_boot

    boot_result.num_boot = num_boot;
    boot_result.bootsamp = reorder;

    if ~isbehav

        sal_sum2 = (sal_sum.^2) / (num_boot);
        dsal_sum2 = (dsal_sum.^2) / (num_boot);

        boot_result.original_sal = original_sal;
        boot_result.brain_se = sqrt((sal_sq - sal_sum2)/(num_boot-1));

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

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

        %  now compare the original unstandardized saliances with the
        %  bootstrap salances
        %
        boot_result.compare = original_sal ./ boot_result.brain_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

        boot_result.orig_designlv = orig_behavlv;
        boot_result.design_se = sqrt((dsal_sq - dsal_sum2)/(num_boot-1));

    else

        sal_sum2 = (sal_sum.^2) / (num_boot+1);
        dsal_sum2 = (dsal_sum.^2) / (num_boot+1);

        %  compute standard errors - standard deviation of bootstrap sample
        %  since original sample is part of bootstrap, divide by number of
        %  bootstrap iterations rather than number of bootstraps minus 1
        %
        %  add ceiling to calculations to prevent the following operations
        %  from producing negative/complex numbers
        %
        brain_se = sqrt((ceil(sal_sq)-ceil(sal_sum2))/(num_boot));
        behav_se = sqrt((ceil(dsal_sq)-ceil(dsal_sum2))/(num_boot));

        %  check for zero standard errors - replace with ones
        %
        test_zeros=find(brain_se<=0);
        if ~isempty(test_zeros);
            brain_se(test_zeros)=1;
        end

        compare = original_sal ./ brain_se;
        compare_behavlv = orig_behavlv ./ behav_se;

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

        ul=Clim;
        ll=100-Clim;

        % e.g. 0.05 >> 0.025 for upper & lower tails, two-tailed
        %
        ClimNi = 0.5*(1-(Clim*0.01));

        %  loop to calculate upper and lower CI limits
        %
        for r=1:r1
            for c=1:c1
                ulcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ul);
                llcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ll);
                prop(r,c)=length( find(distrib(r,c,2:num_boot+1) ...
                                       <= orig_corr(r,c)) ) / num_boot;
                if prop(r,c)==1 |prop(r,c)==0	% can't calculate the cumulative_gaussian_inv
                    llcorr_adj(r,c)=NaN;
                    ulcorr_adj(r,c)=NaN;
                else

                    % adjusted confidence intervals - in case the
                    % bootstrap samples are extremely skewed

                    % norm inverse to start to adjust conf int
                    %
                    ni=cumulative_gaussian_inv(prop(r,c));

                    % 1st part of recalc the lower conf interval,
                    % this evaluates to +1.96 for 95%CI
                    %
                    uli=(2*ni) + cumulative_gaussian_inv(1-ClimNi);

                    % 1st part of recalc the upper conf interval
                    % e.g -1.96 for 95%CI
                    %
                    lli=(2*ni) + cumulative_gaussian_inv(ClimNi);

                    ncdf_lli=cumulative_gaussian(lli)*100;	% percentile for lower bounds
                    ncdf_uli=cumulative_gaussian(uli)*100;	% percentile for upper bounds

                    % new percentile
                    %
                    llcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ...
			ncdf_lli));
                    ulcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ...
			ncdf_uli));

                end	% if
            end		% for c
        end		% for r

        boot_result.orig_corr = orig_corr;
        boot_result.ulcorr = ulcorr;
        boot_result.llcorr = llcorr;
        boot_result.ulcorr_adj = ulcorr_adj;
        boot_result.llcorr_adj = llcorr_adj;
        boot_result.prop = prop;
        boot_result.distrib = distrib;
        boot_result.orig_brainlv = original_sal;
        boot_result.brain_se = brain_se;
        boot_result.compare = compare;
        boot_result.compare_behavlv = compare_behavlv;
        boot_result.badbeh = badbeh;
        boot_result.countnewboot = countnewboot;

    end		% enf if ~isbehav

    return;					% rri_analysis_boot


