
% should be started from the right directory
% ------------------------------------------


[foo,ls_dirs]=unix('ls -d *_0.2');
idx=findstr(ls_dirs,'_0.2')+3;
n_files=length(idx);

roi_dirs=[];
roi_dirs{1}=ls_dirs(1:idx(1));

for i = 2:n_files
  roi_dirs{i}=ls_dirs(idx(i-1)+2:idx(i));
end



file=['No01.mat'; 'No03.mat'; 'No05.mat'; 'No06.mat'; 'No07.mat'; ...
      'No08.mat'; 'No10.mat'; 'No11.mat'; 'No12.mat'; 'No13.mat'];

% should already have 'name' and 'roi';
load('/vol/disk87/kalina/rel_pr/ROIs/Roi.mat');

% should already have 'binned', 'rt_binned'
load('/vol/disk87/kalina/rel_pr/rt_bin_36.mat');

% should already have 'binned_inacc','rt_binned_inacc'
%		      'idx_bin_acc', 'idx_bin_inacc';
load('/vol/disk87/kalina/rel_pr/rt_bin_ai_12.mat');

% should already have 'rels'
load('/vol/disk87/kalina/rel_pr/rav_BehStats.mat');

% should already have 'SGlobal{1:10}' and 'SSessGM{1:10}'
load('/vol/disk87/kalina/rel_pr/global.mat');



%===============================================================
% START CYCLING THROUGH ROI DIRECTORIES
for roi_file = 1:n_files

  cd(char(roi_dirs{roi_file}));



% get the name of the (only) image files in this directory
% (there should be only one .img file)
% --------------------------------------------------------
roi_name = ls('*img');
roi_name = roi_name(1:end-5);
roi_name(find(roi_name=='_'))='-';

roi_idx       = length(name)+1;
name{roi_idx} = roi_name;

for i = 1:10
    load(file(i,:));   
    SubjCond{i}=Cond;
    clear Cond;
end


   % extract data from individual subject's mat files
   % -----------------------------------------------
   for c=1:3

       roi(roi_idx).Y{c}=[];
       for s=1:10
         roi(roi_idx).Y{c} = [roi(roi_idx).Y{c}; SubjCond{s}{c}.rawY'];
       end

       roi(roi_idx).Yzsc{c}=[];
       for s=1:10
         roi(roi_idx).Yzsc{c} = [roi(roi_idx).Yzsc{c}; SubjCond{s}{c}.Yzsc'];
       end


       roi(roi_idx).Ytrm{c}=[];
       for s=1:10
         roi(roi_idx).Ytrm{c} = [roi(roi_idx).Ytrm{c}; SubjCond{s}{c}.Ytrm'];
       end


       roi(roi_idx).Ydz{c}=[];
       for s=1:10
         roi(roi_idx).Ydz{c} = [roi(roi_idx).Ydz{c}; SubjCond{s}{c}.Vdz'];
       end

       roi(roi_idx).Ydztr{c}=[];
       for s=1:10
         roi(roi_idx).Ydztr{c} = [roi(roi_idx).Ydztr{c}; SubjCond{s}{c}.Vdztr'];

       end

       % Ydztr_bin
       % ------------------------------------
       roi(roi_idx).Ydztr_bin{c}=roi(roi_idx).Ydztr{c}(find(binned{c}(:)),:);
  
       roi(roi_idx).Ydztr_bin_tr{c} = roi(roi_idx).Ydztr_bin{c};
      
       roi(roi_idx).Ydztr_bin_tr{c}( ...
                 abs(roi(roi_idx).Ydztr_bin_tr{c})>2);

    end


       % Ydztr_binai for 2-rel
       % ------------------------------------
       c=3;
       roi(roi_idx).Ydztr_binai{c}= ...
             roi(roi_idx).Ydztr{c}(find(binned_inacc{c}(:)),:);
       
     

       % reconstruct the raw time course, by knowing 'rels'
       % 'rels' in this case is the condition in index
       % -------------------------------------------------
       for s=1:10
               Subj{s}.Yraw=[];
               Subj{s}.Ydz =[];

               tmp=ones(1,3);% the trial index for this condition
                             % ranging from 1 to 24

         for i=1:72
  	       Subj{s}.Yraw=[Subj{s}.Yraw; ...
                             SubjCond{s}{rels(i)}.rawY(:,tmp(rels(i)))];
               Subj{s}.Ydz =[Subj{s}.Ydz; ...
                             SubjCond{s}{rels(i)}.Vdz(:,tmp(rels(i)))];
          
               tmp(rels(i))=tmp(rels(i))+1; % update trial index
         end
       end



%%A way to apply global scaling (commented out)
%%for s=1:10
%%   Subj{s}.Yraw =  Subj{s}.Yraw ./ SGlobal{s};
%%end



                %------ use spm's spm_make_filter-------------- 
%                filterLF.Choice='specify';
%                filterLF.Param=64;   
%                filterHF.Choice='none';
%                filterHF.Param=[]; %FWHM (secs)
%                % from spm_results_ui.m, needs some transformations
%                filterHF.Param=filterHF.Param/sqrt(8*log(2));

%                k=1080; % # scans
%                RT=2;   % repetition time in sec

%                [K] = spm_make_filter(k,RT,filterHF,filterLF);

       
       % start pre-processing from the raw data
       % --------------------------------------
 
       for s=1:10

               % FIRST ----------------------------------------
               % remove scanner spikes (mainly from insufficient
               % nextras) by replacing abs(zscore)>2 with session
               % specific trimmed mean
               
               out1 = find(abs(zscore(Subj{s}.Yraw(1:360)))>2);
               out2 = find(abs(zscore(Subj{s}.Yraw(361:720)))>2);
               out3 = find(abs(zscore(Subj{s}.Yraw(720:1080)))>2);

               trM1 = trimmean(Subj{s}.Yraw(1:360),10);
               trM2 = trimmean(Subj{s}.Yraw(361:720),10);
               trM3 = trimmean(Subj{s}.Yraw(721:1080),10);

               Subj{s}.Yraw_tr           = Subj{s}.Yraw;
               Subj{s}.Yraw_tr(out1)     = trM1;
               Subj{s}.Yraw_tr(360+out2) = trM2;
               Subj{s}.Yraw_tr(720+out3) = trM3;
               
               
               % FIRST.AND.A.HALF ----------------------------------
               % apply session specific grand-mean scaling, dividing
               % by the grand mean for the entire brain for this
               % session
                 Subj{s}.Yraw_tr(1:360)   = ...
                 Subj{s}.Yraw_tr(1:360)*100   / SSessGM{s}(1);

                 Subj{s}.Yraw_tr(361:720) = ...
                 Subj{s}.Yraw_tr(361:720)*100 / SSessGM{s}(2);

                 Subj{s}.Yraw_tr(721:1080)= ...
                 Subj{s}.Yraw_tr(721:1080)*100/ SSessGM{s}(3);


               % SECOND ---------------------------------------
               % detrend by removing linear trend on a session
               % by session basis
               Subj{s}.Yraw_tr_dtr =[...
                  detrend(Subj{s}.Yraw_tr(1:360)); ...
                  detrend(Subj{s}.Yraw_tr(361:720)); ...
                  detrend(Subj{s}.Yraw_tr(721:1080))];
 
%                % Or, alternatively, detrend linearly and non-
%                % linearly, using spm_detrend and a
%                % with a 4th order polynomial
%              Subj{s}.Yraw_tr_dtr =[...
%                  spm_detrend(Subj{s}.Yraw_tr(1:360),1); ...
%                  spm_detrend(Subj{s}.Yraw_tr(361:720),1); ...
%                  spm_detrend(Subj{s}.Yraw_tr(721:1080),1)];
%                % Tried that, but it gives the exact same resuls
%                % as the linear detrending above, implying
%                % perhaps that there are no non-linear trends
%                % within session. Also, any detrending should
%                % always be done within session, or otherwise
%                % anomalous results seem to obtain.
                

               % THIRD ---------------------------------------
               % bandpass filter (try low along, and high alone
               % too). The cutoffs are ( 1/16 Hz - 1/64 Hz ),
               % with Nyquist = 1/4 Hz. Butterworth filter.
                  
	       % low-pass 
               % ---------
               [b,a] = butter(1,(1/20)/(1/4));

%                b=spm_fir(2,[(1/64)/(1/4)],[1 0]);
%                a=1;
               Subj{s}.Yraw_tr_dtr_lp = filter(b,a,Subj{s}.Yraw_tr_dtr);
                  
	       % high-pass 
               % ---------
               [b,a] = butter(1,(1/64)/(1/4),'high');
          
%                b=spm_fir(2,[(1/16)/(1/4)],[0 1]);
%                a=1;               
               Subj{s}.Yraw_tr_dtr_hp = filter(b,a,Subj{s}.Yraw_tr_dtr);

               % tried a couple of band-pass IIR filters, but it's 
               % probably better to use a FIR filter
               % --------------------------------------
% butterworth  
[b,a] = butter(1,[(1/64)/(1/4) , (1/20)/(1/4)]);

% elliptic     [b,a] = ellip(1,1,5,[(1/64)/(1/4) , (1/16)/(1/4)]);
% chebyshev II [b,a] = cheby2(3,20,[(1/64)/(1/4) , (1/16)/(1/4)]);

% FIR
% least square linear-phase FIR filter (using spm_fir.m);
%                b=spm_fir(2,[(1/64)/(1/4) , (1/16)/(1/4)],[0 1 0]);
%                a=1;
                Subj{s}.Yraw_tr_dtr_bp = filter(b,a,Subj{s}.Yraw_tr_dtr);


               %------ use spm's spm_make_filter-------------- 
%                Subj{s}.Yraw_tr_dtr_bp = ...
%                     spm_filter('apply',K,Subj{s}.Yraw_tr_dtr);
%               % Tried that, but spm's filter I think is IIR too
%               % JBP seems to have something planned in spm_make_filter
%               % on FIR filtering, but it doesn't really work for now


         end
         % end preprocessing of the raw data



         % Inset the new stuff into the SubjCond structure
         % ---------------------------------------------------
         for s=1:10
         for c=1:3

 	 SubjCond{s}{c}.Yraw_tr =        zeros(15,24);
	 SubjCond{s}{c}.Yraw_tr_dtr =    zeros(15,24);  
         SubjCond{s}{c}.Yraw_tr_dtr_lp = zeros(15,24);    
         SubjCond{s}{c}.Yraw_tr_dtr_hp = zeros(15,24);
         SubjCond{s}{c}.Yraw_tr_dtr_bp = zeros(15,24);

	 SubjCond{s}{c}.Yraw_tr_avg =        zeros(15,1);   
	 SubjCond{s}{c}.Yraw_tr_dtr_avg =    zeros(15,1);   
         SubjCond{s}{c}.Yraw_tr_dtr_lp_avg = zeros(15,1);      
         SubjCond{s}{c}.Yraw_tr_dtr_hp_avg = zeros(15,1);  
         SubjCond{s}{c}.Yraw_tr_dtr_bp_avg = zeros(15,1); 

         

           for i = 1:15    % (num_scans in PST)

             SubjCond{s}{c}.Yraw_tr(i,:) = ...
                    Subj{s}.Yraw_tr(SubjCond{s}{c}.pst_idx == i);
    
	     SubjCond{s}{c}.Yraw_tr_dtr(i,:) =  ...  
                    Subj{s}.Yraw_tr_dtr(SubjCond{s}{c}.pst_idx == i);
    
	     SubjCond{s}{c}.Yraw_tr_dtr_lp(i,:) = ...  
                    Subj{s}.Yraw_tr_dtr_lp(SubjCond{s}{c}.pst_idx == i);

    	     SubjCond{s}{c}.Yraw_tr_dtr_hp(i,:) = ...
                    Subj{s}.Yraw_tr_dtr_hp(SubjCond{s}{c}.pst_idx == i);

    	     SubjCond{s}{c}.Yraw_tr_dtr_bp(i,:) = ...
                    Subj{s}.Yraw_tr_dtr_bp(SubjCond{s}{c}.pst_idx == i);

    
             SubjCond{s}{c}.Yraw_tr_avg(i) = ...
               mean(Subj{s}.Yraw_tr(SubjCond{s}{c}.pst_idx == i));
    
	     SubjCond{s}{c}.Yraw_tr_dtr_avg(i) =  ...  
               mean(Subj{s}.Yraw_tr_dtr(SubjCond{s}{c}.pst_idx == i));
    
	     SubjCond{s}{c}.Yraw_tr_dtr_lp_avg(i) = ...  
               mean(Subj{s}.Yraw_tr_dtr_lp(SubjCond{s}{c}.pst_idx == i));

    	     SubjCond{s}{c}.Yraw_tr_dtr_hp_avg(i) = ...
               mean(Subj{s}.Yraw_tr_dtr_hp(SubjCond{s}{c}.pst_idx == i));

    	     SubjCond{s}{c}.Yraw_tr_dtr_bp_avg(i) = ...
               mean(Subj{s}.Yraw_tr_dtr_bp(SubjCond{s}{c}.pst_idx == i));

           end
         end
         end


     % Insert the new stuff into the roi structure
     % -------------------------------------------
     for c=1:3

         roi(roi_idx).Yraw_tr{c} =  [];
	 roi(roi_idx).Yraw_tr_dtr{c} = [];    
         roi(roi_idx).Yraw_tr_dtr_lp{c} = [];    
         roi(roi_idx).Yraw_tr_dtr_hp{c} = [];
         roi(roi_idx).Yraw_tr_dtr_bp{c} = [];

         for s=1:10

	     roi(roi_idx).Yraw_tr{c} =  [roi(roi_idx).Yraw_tr{c}; ...
                                       SubjCond{s}{c}.Yraw_tr'];
    
	     roi(roi_idx).Yraw_tr_dtr{c} = [roi(roi_idx).Yraw_tr_dtr{c};
                                          SubjCond{s}{c}.Yraw_tr_dtr'];
    
	     roi(roi_idx).Yraw_tr_dtr_lp{c} = [roi(roi_idx).Yraw_tr_dtr_lp{c};
                                             SubjCond{s}{c}.Yraw_tr_dtr_lp'];
    
	     roi(roi_idx).Yraw_tr_dtr_hp{c} = [roi(roi_idx).Yraw_tr_dtr_hp{c};
                                             SubjCond{s}{c}.Yraw_tr_dtr_hp'];
    
             roi(roi_idx).Yraw_tr_dtr_bp{c} = [roi(roi_idx).Yraw_tr_dtr_bp{c};
                                             SubjCond{s}{c}.Yraw_tr_dtr_bp'];

         end 

     end




     % BIN the added stuff to equate for RT across conditions
     % -------------------------------------------------------
     for c=1:3
       roi(roi_idx).Yraw_tr_dtr_bp_bin{c} = ...
       roi(roi_idx).Yraw_tr_dtr_bp{c}(find(binned{c}(:)),:);        
     end

       % acc vs. inacc for 2-rel
       c=3;
       roi(roi_idx).Yraw_tr_dtr_bp_binai{c} = ...
       roi(roi_idx).Yraw_tr_dtr_bp{c}(find(binned_inacc{c}(:)),:);



        % Apply global scaling to see what would happen;
         % ----------------------------------------------
	   for s=1:10
   	     Subj{s}.Yraw_gs =  Subj{s}.Yraw ./ SGlobal{s};
           end


           for s=1:10
           for c=1:3

             SubjCond{s}{c}.Yraw_gs =        zeros(15,24);
             SubjCond{s}{c}.Yraw_gs_avg =    zeros(15,1);

            for i = 1:15    % (num_scans in PST)
             SubjCond{s}{c}.Yraw_gs(i,:) = ...
                    Subj{s}.Yraw_gs(SubjCond{s}{c}.pst_idx == i);
             SubjCond{s}{c}.Yraw_gs_avg(i) = ...
               mean(Subj{s}.Yraw_gs(SubjCond{s}{c}.pst_idx == i));
            end
     
           end
           end
      

           for c=1:3             
             roi(roi_idx).Yraw_gs{c} = [];
                 for s=1:10
	         roi(roi_idx).Yraw_gs{c} =  [roi(roi_idx).Yraw_gs{c}; ...
                                           SubjCond{s}{c}.Yraw_gs'];
                 end
           end

          % end global scaling application
          % ------------------------------

% save Roi roi name rt acc binned rt_binned binned_inacc rt_binned_inacc;    

    leg = ['0rel';'1rel';'2rel'];
    pst = 0:2:28;
    COL   = {'g','b','r'};
   

%    subplot(3,4,1);
%    title(name{roi_idx});
%    hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Y{c},1),char(COL(c)),'LineWidth',1)
%          legend(leg);
%         end

%     subplot(3,4,2);
%     title('trimmed')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Ytrm{c},1),char(COL(c)),'LineWidth',1)
%         end


%     subplot(3,4,3);
%     axis([0 28 -0.8 1]);
%     title('zsc trm')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Yzsc{c},1),char(COL(c)),'LineWidth',1)
%         end


%     subplot(3,4,4);
%     axis([0 28 -0.8 1]);
%     title('BIN trimmed')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Ydztr_bin_tr{c},1),char(COL(c)),'LineWidth',1)
%         end

     
%    subplot(3,4,5);
%     axis([0 28 -0.8 1]);
%     title('detr filt zsc trm')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Ydztr{c},1),char(COL(c)),'LineWidth',1)
%         end


%     subplot(3,4,6);
%     axis([0 28 -0.8 1]);
%     title('BIN detr filt zsc trm')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Ydztr_bin{c},1),char(COL(c)),'LineWidth',1)
%         end

%     subplot(2,3,1);
%     title(name{roi_idx});
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Yraw_tr{c},1),char(COL(c)),'LineWidth',)
%         end

%     subplot(2,3,2);
%     axis([0 28 -0.8 1]);
%     title('raw tr detr')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Yraw_tr_dtr{c},1),char(COL(c)),'LineWidth',1)
%         end


%     subplot(2,3,3);
%     axis([0 28 -0.8 1]);
%     title('raw tr detr LP')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Yraw_tr_dtr_lp{c},1),char(COL(c)),'LineWidth',1)
%         end

%     subplot(2,3,4);
%     axis([0 28 -0.8 1]);
%     title('raw tr detr HP')
%     hold on;
%         for c = 1:3
%           plot(pst,mean(roi(roi_idx).Yraw_tr_dtr_hp{c},1),char(COL(c)),'LineWidth',1)
%         end


% END CYCLING THROUGH ROI DIRS
cd ../;
end



% Plot all for all rois
% ---------------------
h1=figure; hold on;
for roi_file = 1:n_files
     
     roi_idx=roi_file+1;
    
     tit= char(roi_dirs{roi_file});
     tit([findstr(tit,'_')])='-';
     tit=[tit,' ALL'];
     subplot(2,3,roi_file);
     axis([0 28 -0.8 0.8]);
     title(tit);
     hold on;
 
         for c = 1:3
           mean_pst=trimmean(roi(roi_idx).Yraw_tr_dtr_bp{c},10);
           plot(pst,mean_pst,char(COL(c)),'LineWidth',2);

           peak_idx=find(mean_pst==max(mean_pst));
             if peak_idx<4 | peak_idx>10
                 peak_idx=find(mean_pst==min(mean_pst));
             end
       

           % SEM computation
           s=1:24:217;
           for i=1:10
             SubjMean(i)=  mean(roi(roi_idx).Yraw_tr_dtr_bp{c}...
                                           (s(i):s(i)+23,peak_idx)); 
           end

           SEM  =std(SubjMean)/...
                   sqrt(10);

           errorbar_k(pst(peak_idx),mean_pst(peak_idx),...
                      SEM,char(COL(c)));
         end
end
%------ 


%--------plot binned for all rois
h2=figure; hold on;
for roi_file = 1:n_files
     
     roi_idx=roi_file+1;

     tit= char(roi_dirs{roi_file});
     tit([findstr(tit,'_')])='-';
     tit=[tit,' BIN'];
     subplot(2,3,roi_file);
     axis([0 28 -0.8 0.8]);
     title(tit);
     hold on;
 
         for c = 1:3
           mean_pst=trimmean(roi(roi_idx).Yraw_tr_dtr_bp_bin{c},10);
           plot(pst,mean_pst,char(COL(c)),'LineWidth',2);

           peak_idx=find(mean_pst==max(mean_pst));
             if peak_idx<4 | peak_idx>10
                 peak_idx=find(mean_pst==min(mean_pst));
             end
       

           SEM  =std(roi(roi_idx).Yraw_tr_dtr_bp_bin{c}(:,peak_idx))/...
                   sqrt(36);
%           errorbar_k(pst(peak_idx),mean_pst(peak_idx),...
%                      SEM,char(COL(c)));
         end
end
%---


%--------plot acc vs. inacc binned for all rois
h3=figure; hold on;
for roi_file = 1:n_files
     
     roi_idx=roi_file+1;

     tit= char(roi_dirs{roi_file});
     tit([findstr(tit,'_')])='-';
     tit=[tit,' BIN AI'];
     subplot(2,3,roi_file);
     axis([0 28 -0.8 0.8]);
     title(tit);
     hold on;
 
         c=3;

         tmp_plot{1}=roi(roi_idx).Yraw_tr_dtr_bp_binai{c}(...
                                     idx_bin_acc,:);
         tmp_plot{2}=roi(roi_idx).Yraw_tr_dtr_bp_binai{c}(...
                                     idx_bin_inacc,:);
 
         % use the trimmed mean, 5% in each tail (z=2)
         COLai = {'r--','r-.'};
         for a=1:2
           mean_pst=trimmean(tmp_plot{a},10);
           plot(pst,mean_pst,char(COLai(a)),'LineWidth',3);

           peak_idx=find(mean_pst==max(mean_pst));
             if peak_idx<4 | peak_idx>10
                 peak_idx=find(mean_pst==min(mean_pst));
             end
       

           SEM  =std(tmp_plot{a}(:,peak_idx))/...
                   sqrt(12);
%           errorbar_k(pst(peak_idx),mean_pst(peak_idx),...
%                      SEM,char(COLai(a)));

         end  
end
%----


















