% get current directry
cwd=pwd;

ROI_file=spm_get(Inf, '.mat', 'Select ROI.mat file(s) to plot');

if exist('plot_info.mat') == 2
  load plot_info;
end

%--------------------------------------
%  Collect details needed for averaging
%--------------------------------------

     %---- define the peristimulus time ----
      if exist('time')==0;

       num_scans=input('How many scans per trial: ');
       num_trials=input('How many trials in the experiment: ');
       pst=[1:num_scans];
      
       time=[];
        for i=1:num_trials     
         time=[time pst];  
        end

        %-Quick check
	fprintf('\nThere should have been ');
        fprintf('%8d',size(time,2));
        fprintf(' scans altogether in the experiment.\n')
        fprintf('\n');
      
      end 

      %-------- define condition types -------
      if exist('trial_type')==0;
     
         %-ask for names of the different trial types
         num_types=input('How many trial types (conditions): ');

             for i=1:num_types
              fprintf('Type'); fprintf('%2d',i);
              trial_type(i,:)=input([' - please enter name: '],'s');
             end
      end


      %-------- define trial types -----------
      if exist('vec_types')==0;
         
         fprintf('\n\nPlease, enter (the name of a matlab variable ');
         fprintf('representing) \na vector of length ');
         fprintf('%4d',num_trials);
         fprintf(' (or the number of trials). \n\nThe vector is needed in ');
         fprintf('order to define each trial as to the \ncondition ');
         fprintf('it belongs to. \nEnter 1 for '); fprintf(trial_type(1,:)); 
         fprintf(', 2 for '); fprintf(trial_type(2,:)); 
         fprintf(' and so forth.\n\n');
   
         vec_types=input('Vector defining the types: ');

      end

      %---------------     
      % replicate the vec_types to form an index vector, identifying
      % each scan as to which scondition it belongs to
      %-------------------------------
      type_idx = reshape(repmat(vec_types,num_scans,1),1,size(time,2));

      
      % Save plotting information in plot_info.mat
      cd(cwd);      
      save plot_info num* pst time trial_type vec_types type_idx;


%======================================================
% start cycling through all mat_files
%=====================================================

  for mat_file = 1:size(ROI_file,1)

    load(ROI_file(mat_file,:));

      %-----------------------------------------------------
      %-average Yc and y across time
      %-----------------------------------------------------

      for scan=1:num_scans                    
         for cond=1:num_types

            t_avgROI.Y(scan,cond)=...
                    mean(avgROI.Y  (find((time==scan)&(type_idx==cond)))); 

            t_avgROI.Yc(scan,cond)=...
                    mean(avgROI.Yc (find((time==scan)&(type_idx==cond)))); 

            t_avgROI.y(scan,cond)=...
                    mean(avgROI.y  (find((time==scan)&(type_idx==cond)))); 

            t_avgROI.Yng(scan,cond)=...
                    mean(avgROI.Yng(find((time==scan)&(type_idx==cond))));
                   
         end
      end



      %----------------------------------------------------- 
      %-plot fitted response, and adjusted data +/- SEM
      %----------------------------------------------------- 


      COL   = ['r','b','g','c','y','m'];

     figure
     hold on
     axis([0 max(pst) floor(min(min(t_avgROI.y)))  ceil(max(max(t_avgROI.y)))])

      for i=1:num_types
            plot(pst,t_avgROI.y(:,i),[COL(i) '-'],'LineWidth',3)
            legend(trial_type)
      end

      for i=1:num_types
            plot(pst,t_avgROI.Yc(:,i),[COL(i) '--'],'LineWidth',1)
      end

%      errorbar(pst,t_avgROI.y(:,1),t_avgROI.SEM_y(:,1),'g')
%      errorbar(pst,t_avgROI.y(:,2),t_avgROI.SEM_y(:,2),'b')
%      errorbar(pst,t_avgROI.y(:,3),t_avgROI.SEM_y(:,3),'r')

      title(['Subject averaged Fit and Adj data for ' ROIdat{1}.name]);

      print -dpsc -append -r300 figures;

%    figure
%    hold on
%    axis([0 max(pst) floor(min(min(t_avgROI.Y)))  ceil(max(max(t_avgROI.Y)))])
      
%     for i=1:num_types
%      plot(pst,t_avgROI.Y(:,i),[COL(i) '-'],'LineWidth',2)
%      legend(trial_type)
%     end

%     title(['Subject averaged Raw (scaled) data for ' ROIdat{1}.name]);


%   figure
%   hold on
%   axis([0 max(pst) floor(min(min(t_avgROI.Yng)))...
%                     ceil(max(max(t_avgROI.Yng)))])

%     for i=1:num_types
%      plot(pst,t_avgROI.Yng(:,i),[COL(i) '-'],'LineWidth',3)
%      legend(trial_type)
%     end

%     title(['Subject averaged Raw (non-scaled) data for ' ROIdat{1}.name]);

%update the BlockData{} matrix
EventData{mat_file}      =t_avgROI;
EventData{mat_file}.name =ROIdat{1}.name;

%clear up newly created variables before cycling through
%the next mat_file
clear ROIdat avgROI t_avgROI

%=======================================================
% end cycling through mat_files
%=======================================================
end





%%%%%%




