%PLS_ANALYSIS Run PLS analysis on one or more datamats. Datamat list cell
%	array, number of subject list array, and number of condition must
%	be input to the program.
%
%	If this is a non-rotate task PLS analysis or a behavior analysis, 
%	contrast text file or behavior text file must also be available. 
%
%	For contrast text file, columns stand for designs, number of rows
%	for each group is equal to number of conditions. Since all the 
%	groups are stacked vertically, total number of rows in contrast 
%	text file should equal to number of groups multiplied by number 
%	of conditions. 
%
%	For behavior text file, columns stand for measurements, number of
%	rows for each group is the same as the datamat for that group.
%	All the behavior group should also be stacked together vertically
%	with a total number equal to the sum of the number of rows in
%	datamat list array.
%
%  Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond);
%
%  Here:
%
%	First input parameter 'datamat_lst' is a datamat cell array. Each
%	cell stands for one datamat, which is also referred as a group. 
%	All datamats must be in the form of "subject in condition".
%
%	Input 'num_subj_lst' is an array contains the number of subjects
%	in each group.
%
%	Input 'num_cond' is the number of conditions in datamat_lst.
%
%	Output 'result' is a struct contains all necessary items that are
%	generated based on different scenarios. scenario question will be
%	asked at the beginning of the analysis.
%
%	The 'result' struct could have items like:
%
%       u:			brainlv or salience
%	s:			singular value
%	v:			designlv or behavlv
%	usc:			brainscores or scalpscores
%	vsc:			designscores or behavscores
%	datamatcorrs_lst:	correlation of behavior data with datamat,
%				only available in behavior PLS.
%	lvcorrs:		correlation of behavior data with usc,
%				only available in behavior PLS.
%	origpost:		posthoc result, only available in behavior
%				PLS, when posthoc file is provided.
%	perm_result:		struct containing permutation result
%		num_perm:	number of permutation
%		sp:		permuted singular value greater than observed
%		sprob:		sp normalized by num_perm
%	boot_result:		struct containing bootstrap result
%		num_boot:	number of bootstrap
%		bootsamp:	bootstrap reorder sample
%		orig_u:		original salience (orig_u = u * s)
%		[u]:		same as u in non-rotate task PLS analysis
%		orig_v:		original designlv or original behavlv
%		[v]:		same as v in non-rotate task PLS analysis
%		compare_u:	compared salience or compared brain
%		compare_v:	compared designlv or compared behavlv
%		u_se:		standard error of salience or brainlv
%		v_se:		standard error of designlv or behavlv
%
%		following boot_result only available in behavior PLS:
%
%		origcorr:	same as lvcorrs
%		ulcorr:		upper boundary of lvcorrs
%		llcorr:		lower boundary of lvcorrs
%		prop:		lvcorrs probability
%		distrib:	lvcorrs distribution
%		ulcorr_adj:	percentile of lvcorrs distribution with
%				lower boundary of lvcorrs
%		llcorr_adj:	percentile of lvcorrs distribution with
%				lower boundary of lvcorrs
%		badbeh:		display bad behav data that is caused by
%				bad re-order (with 0 standard deviation)
%				which will further cause divided by 0
%		countnewboot:	count the new bootstrap sample that is
%				re-ordered for badbeh

%   Created on 05-JAN-2005 by Jimmy Shen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function result = pls_analysis(datamat_lst, num_subj_lst, k)

   result = [];

   %  check input variable
   %
   if nargin ~= 3
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond);');
      disp(' ');disp(' ');
      return;
   end

   if ~iscell(datamat_lst)
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond);');
      disp(' ');
      disp('datamat_lst should be a datamat cell array. Each cell');
      disp('represents a datamat. All datamats must be in the form');
      disp('of "subject in condition".');
      disp(' ');disp(' ');
      return;
   end

   if ~isnumeric(num_subj_lst)
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond);');
      disp(' ');
      disp('num_subj_lst should be a numeric array. Each element');
      disp('represents the number of subjects in the corresponding');
      disp('datamat.');
      disp(' ');disp(' ');
      return;
   end

   if ~isnumeric(k) | length(k) ~= 1
      disp(' ');
      disp('Usage: result = pls_analysis(datamat_lst, num_subj_lst, num_cond);');
      disp(' ');
      disp('num_cond should be a number represents the number of conditions.');
      disp(' ');disp(' ');
      return;
   end

   %  init
   %
   num_groups = length(datamat_lst);
   total_rows = 0;

   %  total rows in stacked datamat (extract datamat from each group,
   %  and stacked them together)
   %
   for g = 1:num_groups
      total_rows = total_rows + size(datamat_lst{g}, 1);
   end

   %  which PLS analysis do you choose to run
   %
   clc;
   method = input('PLS Option:\n1. Mean-Centering Task PLS\n2. Non-Rotated Task PLS\n3. Behavior PLS\nPlease select: ','s');
   method = str2num(method);
   while isempty(method) | (method ~= 1 & method ~= 2 & method ~= 3)
      disp(' ');
      disp('Please select either 1, 2 or 3 and then strike ''Enter''');
      method = input('PLS Option:\n1. Mean-Centering Task PLS\n2. Non-Rotated Task PLS\n3. Behavior PLS\nPlease select: ','s');
      method = str2num(method);
   end

   %  For test purpose
   %
   if exist('plslog.m','file')
      switch method
      case 1
         plslog('Run Mean-Centering Task PLS Analysis');
      case 2
         plslog('Run Non-Rotated Task PLS Analysis');
      case 3
         plslog('Run Behavior PLS Analysis');
      end
   end

   %  number of permutation
   %
   clc;
   num_perm = input('Please enter number of permutation [0]: ','s');
   num_perm = str2num(num_perm);
   if isempty(num_perm) | num_perm <= 0
      num_perm = 0;
   end;

   %  number of bootstrap
   %
   clc;
   num_boot = input('Please enter number of bootstrap [0]: ','s');
   num_boot = str2num(num_boot);
   if isempty(num_boot) | num_boot <= 0
      num_boot = 0;
   end;

   %  default value
   %
   contrast_file = '';
   behav_file = '';
   posthoc_file = '';
   Clim = 95;

   %  need contrast file for non-rotate task PLS analysis
   %
   if method == 2
      clc;
      done = 0;
      while ~done
         contrast_file = input('Contrast filename with path for Non-rotate PLS: ','s');

         if ~isempty(contrast_file) & exist(contrast_file,'file')
            try
               stacked_designdata = load(contrast_file);
            catch
               disp(' '); disp('Invalid contrast data file.');
               return;
            end

            if size(stacked_designdata,1) ~= num_groups * k
               disp(' '); disp('Wrong number of rows in contrast data file, please try again');
            else
               done = 1;
            end
         else
            disp(' '); disp('File does not exist, please try again');
         end
      end
   end

   %  need behavior file and posthoc file for behavior PLS analysis
   %
   if method == 3
      clc;
      done = 0;
      while ~done
         behav_file = input('Behav filename with path for Behavior PLS: ','s');

         if ~isempty(behav_file) & exist(behav_file,'file')
            try
               stacked_behavdata = load(behav_file);
            catch
               disp(' '); disp('Invalid behavior data file.');
               return;
            end

            if size(stacked_behavdata,1) ~= total_rows
               disp(' '); disp('Wrong number of rows in behavior data file, please try again');
            else
               done = 1;
            end
         else
            disp(' '); disp('File does not exist, please try again');
         end
      end

      %  posthoc text file
      %
      clc;
      posthoc_file = input('Posthoc filename with path for Behavior PLS []: ','s');
      while ~isempty(posthoc_file) & ~exist(posthoc_file,'file')
         disp(' '); disp('Please enter a valid file or strike ''Enter'' if there is no posthoc file');
         posthoc_file = input('Posthoc filename with path for Behavior PLS []: ','s');
      end

      if ~isempty(posthoc_file) & exist(posthoc_file,'file')
         try
            posthoc = load(posthoc_file);

            if size(posthoc,1) ~= size(stacked_behavdata,2)*k*num_groups
               disp(' '); disp('Invalid posthoc row number, posthoc computation will not be performed');
               posthoc_file = '';
               hosthoc = [];
               disp(' '); disp('Press any key to continue');
               pause;
            end
         catch
            disp(' '); disp('Invalid posthoc data file, posthoc computation will not be performed');
            posthoc_file = '';
            hosthoc = [];
            disp(' '); disp('Press any key to continue');
            pause;
         end
      end

      %  confidence level
      %
      clc;
      if num_boot ~= 0
         Clim = input('Please enter a confidence level for Behavior PLS [95]: ','s');
         Clim = str2num(Clim);
         if isempty(Clim) | Clim < 0 | Clim > 100
            Clim = 95;
         end;
      end
   end

   %  init variable for the following loop
   %
   vsc = [];
   datamatcorrs_lst = {};
   stacked_datamat = [];
   stacked_datamatcorrs = [];
   first = 1;			% first row for stacked_behavdata

   %  loop accross the groups, and
   %  calculate datamatcorrs for each group
   %
   for g = 1:num_groups
      clc;
      disp(['Stacking datamat from group: ', num2str(g)]);

      datamat = datamat_lst{g};
      stacked_datamat = [stacked_datamat; datamat];

      n = num_subj_lst(g);

      % compute correlation or covariance
      %
      switch method
      case 1
%         datamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat); % pet & erp style
         meanmat = rri_task_mean(datamat,n);
         datamatcorrs = meanmat-(ones(k,1)*mean(meanmat));
      case 2
         datamatcorrs = rri_task_mean(datamat,n);
      case 3
         last = first + size(datamat,1) - 1;
         behavdata = stacked_behavdata(first:last, :);
         first = last + 1;
         datamatcorrs = rri_corr_maps(behavdata, datamat, n, k);
         datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}];
      end

      stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs];
   end

   %  save datamatcorrs_lst, if behav PLS
   %
   if method == 3
      result.datamatcorrs_lst = datamatcorrs_lst;
   end

   clc;
   disp('Calculating LVs ...');

   if method == 2	% different computation for non-rotate task PLS

      crossblock = rri_normalize(stacked_designdata)'*stacked_datamatcorrs;
      u = crossblock';
      s = sqrt(sum(crossblock.^2, 2));
      v = stacked_designdata;
      normalized_u = rri_normalize(u);
      lvintercorrs = normalized_u'*normalized_u;
      result.lvintercorrs = lvintercorrs;

   else			% mean-centering and behavior PLS

      %  Singular Value Decomposition, observed
      %
      [r c] = size(stacked_datamatcorrs);
      if r <= c
         [u,s,v] = svd(stacked_datamatcorrs',0);
      else
         [u,s,v] = svd(stacked_datamatcorrs,0);
      end

      s = diag(s);

      if ~isempty(posthoc_file)
         origpost = rri_xcor(posthoc,v);
         porigpost = zeros(size(origpost));
         result.origpost = origpost;
      end

   end

   %  save u,s,v for all situations
   %
   result.u = u;
   result.s = s;
   result.v = v;

   clc;
   disp('Calculating Scores ...');

   switch method
   case {1, 2}
      usc = stacked_datamat * u;
      num_col = size(v, 2);

      % expand the num_subj for each row (cond)
      % did the samething as testvec
      %
      for g = 1:num_groups
         n = num_subj_lst(g);

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

         vsc = [vsc; tmp];			% stack by groups
      end
   case 3
      [usc, vsc, lvcorrs] = rri_get_behavscores(stacked_datamat, ...
		stacked_behavdata, u, v, k, num_subj_lst);

      result.lvcorrs = lvcorrs;

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

   %  save Scores for all situations
   %
   result.usc = usc;
   result.vsc = vsc;

   %  save Input for all situations
   %
   result.datamat_lst = datamat_lst;
   result.num_subj_lst = num_subj_lst;
   result.num_conditions = k;

   if num_perm > 0
      sp = zeros(size(s));		% permcount

      if method == 3
         for p = 1:num_perm
            reorder(:,p) = [randperm(size(stacked_datamat,1))'];
         end
      else
         reorder = rri_perm_order(num_subj_lst, k, num_perm);
      end

      for p = 1:num_perm
         clc;
         disp(['Working on Permutation: ',num2str(p),' out of ',num2str(num_perm)]);

         if method == 3
            behav_p = stacked_behavdata(reorder(:,p),:);
         else
            data_p = stacked_datamat(reorder(:,p),:);
         end

         stacked_data_p = [];

         for g=1:num_groups
            n = num_subj_lst(g);
            span = sum(num_subj_lst(1:g-1)) * k;

            switch method
            case 1
               if num_groups == 1
                  meanmat = rri_task_mean(data_p, n);
                  data = meanmat - (ones(k,1)*mean(meanmat));
               else
                  meanmat = rri_task_mean(data_p(1+span:n*k+span,:), n);
                  data = meanmat - (ones(k,1)*mean(meanmat));
               end
            case 2
               if num_groups == 1
                  data = rri_task_mean(data_p, n);
               else
                  data = rri_task_mean(data_p(1+span:n*k+span,:), n);
               end
            case 3
		% Check for upcoming NaN and re-sample if necessary.
		% this only happened on behavior analysis, because the
		% 'xcor' inside of 'rri_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.
		%
                min1 = min(std(behav_p(1+span:n*k+span,:)));
                while (min1 == 0)
                    reorder(:,p) = [randperm(size(stacked_datamat,1))'];
                    behav_p = stacked_behavdata(reorder(:,p),:);
                    min1 = min(std(behav_p(1+span:n*k+span,:)));
                end

		% Notice here that stacked_datamat is used, instead of
		% boot_p. This is only for behavpls_perm.
		%
                if num_groups == 1
                   data = rri_corr_maps(behav_p, stacked_datamat, n, k);
                else
                   data = rri_corr_maps(behav_p(1+span:n*k+span,:), ...
				stacked_datamat(1+span:n*k+span,:), n, k);
                end
            end

            stacked_data_p = [stacked_data_p; data];
         end

         if method == 2
            crossblock = rri_normalize(stacked_designdata)'*stacked_data_p;
            sperm = sqrt(sum(crossblock.^2,2));
            sp = sp + (sperm >= s);
         else
            %  Singular Value Decomposition
            %
            [r c] = size(stacked_data_p);
            if r <= c
               [pu, sperm, pv] = svd(stacked_data_p',0);
            else
               [pu, sperm, pv] = svd(stacked_data_p,0);
            end

            %  rotate pv to align with the original v
            %
            rotatemat = rri_bootprocrust(v, pv);

            %  rescale the vectors
            %
            pv = pv * sperm * rotatemat;

            sperm = sqrt(sum(pv.^2));
            sp = sp + (sperm'>=s);

            if ~isempty(posthoc_file)
               tmp = rri_xcor(posthoc, pv);
               porigpost = porigpost + (abs(tmp) >= abs(origpost));
            end
         end		% if method
      end		% for num_perm

      result.perm_result.num_perm = num_perm;
      result.perm_result.permsamp = reorder;
      result.perm_result.sp = sp;
      result.perm_result.sprob = sp ./ num_perm;
   end

   if num_boot > 0
      original_u = u * diag(s);
      original_v = v * diag(s);

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

      %  include original un-resampled order only for mean-centering task PLS
      %
      if method == 1
         [reorder, new_num_boot] = rri_boot_order(num_subj_lst,k,num_boot,1);
      else
         [reorder, new_num_boot] = rri_boot_order(num_subj_lst,k,num_boot);
      end

      if isempty(reorder)
         disp('bootstrap order is not available');
         return;
      end;

      if new_num_boot ~= num_boot
         num_boot = new_num_boot;

         if method == 3
            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

      switch method
      case 1
         u_sq = zeros(size(u));
         u_sum = zeros(size(u));
         v_sq = zeros(size(v));
         v_sum = zeros(size(v));
      case 2
         u_sq = u.^2;
         u_sum = u;
         v_sq = v.^2;
         v_sum = v;
      case 3
         u_sq = original_u.^2;
         u_sum = original_u;
         v_sq = original_v.^2;
         v_sum = original_v;
      end

      for p=1:num_boot
         clc;
         disp(['Working on Bootstrap: ',num2str(p),' out of ',num2str(num_boot)]);

         data_p = stacked_datamat(reorder(:,p),:);
         stacked_data_p = [];

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

            switch method
            case 1
               if num_groups == 1
                  meanmat = rri_task_mean(data_p, n);
                  data = meanmat - (ones(k,1)*mean(meanmat));
               else
                  meanmat = rri_task_mean(data_p(1+span:n*k+span,:), n);
                  data = meanmat - (ones(k,1)*mean(meanmat));
               end
            case 2
               if num_groups == 1
                  data = rri_task_mean(data_p, n);
               else
                  data = rri_task_mean(data_p(1+span:n*k+span,:), n);
               end
            case 3
                if ~is_boot_samples

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

                   % init badbehav array to 0
                   % which is used to record the bad behav data caused by
                   % bad re-order. This variable is for disp only.
                   %
                   badbehav = zeros(k, 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
                           disp(['count_new_boot exceeds num_boot: ', ...
                              num2str(countnewboot)]);
                           breakon=1;
                           break;
                       end

                       reorder(:,p) = rri_boot_order(num_subj_lst, k, 1);

                       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		% switch

            stacked_data_p = [stacked_data_p; data];
         end		% for num_groups

         if method == 2
            crossblock = rri_normalize(stacked_designdata)'*stacked_data_p;
            u_sq = u_sq + (crossblock.^2)';
            u_sum = u_sum + crossblock';
         else
            %  Singular Value Decomposition
            %
            [r c] = size(stacked_data_p);
            if r <= c
               [pu, sboot, pv] = svd(stacked_data_p',0);
            else
               [pu, sboot, pv] = svd(stacked_data_p,0);
            end

            %  rotate pv to align with the original v
            %
            rotatemat = rri_bootprocrust(v, pv);

            %  rescale the vectors
            %
            pu = pu * sboot * rotatemat;
            pv = pv * sboot * rotatemat;

            if method == 3
               [brainsctmp, behavsctmp, bcorr] = ...
			rri_get_behavscores(data_p, behav_p, ...
			rri_normalize(pu), ...
			rri_normalize(pv), k, num_subj_lst);

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

            u_sum = u_sum + pu;
            u_sq = u_sq + pu.^2;

            v_sum = v_sum + pv;
            v_sq = v_sq + pv.^2;
         end
      end		% for num_boot

      result.boot_result.num_boot = num_boot;

      switch method
      case 1

        u_sum2 = (u_sum.^2) / (num_boot);
        v_sum2 = (v_sum.^2) / (num_boot);

        u_se = sqrt((u_sq-u_sum2) / (num_boot-1));
        v_se = sqrt((v_sq-v_sum2) / (num_boot-1));

      case 2

        %  calculate standard error
        %
        u_sum2 = (u_sum.^2) / (num_boot+1);
        v_sum2 = (v_sum.^2) / (num_boot+1);

        u_se = sqrt((u_sq-u_sum2) / num_boot);
        v_se = sqrt((v_sq-v_sum2) / num_boot);

      case 3

        u_sum2 = (u_sum.^2) / (num_boot+1);
        v_sum2 = (v_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
        %
        u_se = sqrt((ceil(u_sq)-ceil(u_sum2))/(num_boot));
        v_se = sqrt((ceil(v_sq)-ceil(v_sum2))/(num_boot));

        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

        result.boot_result.orig_corr = orig_corr;
        result.boot_result.ulcorr = ulcorr;
        result.boot_result.llcorr = llcorr;
        result.boot_result.ulcorr_adj = ulcorr_adj;
        result.boot_result.llcorr_adj = llcorr_adj;
        result.boot_result.prop = prop;
        result.boot_result.distrib = distrib;
        result.boot_result.badbeh = badbeh;
        result.boot_result.countnewboot = countnewboot;

      end

      result.boot_result.bootsamp = reorder;

      % result.boot_result.u_sum2 = u_sum2;
      % result.boot_result.v_sum2 = v_sum2;

      if method == 2
         result.boot_result.u = u;
         result.boot_result.v = v;
      else
         result.boot_result.orig_u = original_u;
         result.boot_result.orig_v = original_v;
      end

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

      test_zeros_v=find(v_se<=0);
      if ~isempty(test_zeros_v);
         v_se(test_zeros_v)=1;
      end

      if method == 2
         compare_u = u ./ u_se;
         compare_v = v ./ v_se;
      else
         compare_u = original_u ./ u_se;
         compare_v = original_v ./ v_se;
      end

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

      if ~isempty(test_zeros_v);
         compare_v(test_zeros_v)=0;
      end

      result.boot_result.compare_u = compare_u;
      result.boot_result.compare_v = compare_v;
      result.boot_result.u_se = u_se;
      result.boot_result.v_se = v_se;

   end