% sig_topo() - Plots scalp topographies of ERPs/ERP efffects that have been
%              tested for significance in a family of tests.  Only works for
%              ERPs/ERP effects derived from averaging time points within 
%              time windows (for tests of every single time point within
%              time windows use gui_erp.m).  Topographies can be visualized
%              in microvolts or t-scores.
%             
% Usage:
%  >> sig_topo(GND_GRP_or_fname,t_test_id,varargin)
%
% Required Inputs:
%  GND_GRP_or_fname - A GND or GRP structure variable or the 
%                      filename of a GND/GRP structure that has
%                      been saved to disk. If you specifiy a filename be
%                      sure to include the file's path, unless the file is
%                      in the current working directory.
%  t_test_id        - [integer] The index # of the set of t-test results 
%                      stored in the GND/GRP variable that you wish to visualize.  
%                      To see what test results are available, look at 
%                      GND.t_tests, GRP.t_tests, or run headinfo(GND).
%
% Optional Inputs:
%  fig_id           - [integer] The index # of the MATLAB figure in which
%                     the diagram will be produced.  Useful for overwriting
%                     old figures. {default: lowest unused index}
%  units            - ['uV' or 't'] If 'uV' topographies will be visualized
%                     in units of microvolts.  If 't', they will be shown 
%                     in t-scores. {default: 'uV'}
%  title_on         - [0 or 1] If 1, the bin # and bin descriptor will be
%                     printed at the top of the figure. {default: 1}
%  one_scale        - [0 or 1] If 1, all topographies will be plot using
%                     the same color scale (useful for comparing effect 
%                     sizes across topographies).  If 0, each topography
%                     will be plot using a color scale tailored to it
%                     (provides maximal resolution of topography
%                     distribution).  {default: 1 if plotting in uV, 0 if 
%                     plotting t-scores}
%  verblevel        - An integer specifiying the amount of information you want
%                     this function to provide about what it is doing during runtime.
%                     Options are:
%                      0 - quiet, only show errors, warnings, and EEGLAB reports
%                      1 - stuff anyone should probably know
%                      2 - stuff you should know the first time you start working
%                          with a data set {default value}
%                      3 - stuff that might help you debug (show all
%                          reports)
%
% Global Variable:
%   VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells 
%               functions how much to report about what they're doing during             
%               runtime) set by the optional function argument 'verblevel'
%
% Notes:
% -The printed/exported figure will have the same dimensions as the figure
% on the display.  Thus you can undo figure clutter by re-sizing it.
%
% Author:
% David Groppe
% Kutaslab, 3/2010

%%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
% 5/6/2010-Made compatible with FDR control

%%%%%%%%%%%%%%%% FUTURE WORK %%%%%%%%%%%%%%%%%
% -

function sig_topo(GND_GRP_or_fname,t_test_id,varargin)

p=inputParser;
p.addRequired('GND_GRP_or_fname',@(x) ischar(x) || isstruct(x));
p.addRequired('t_test_id',@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('units','uV',@ischar);
p.addParamValue('fig_id',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('title_on',1,@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('one_scale',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));

p.parse(GND_GRP_or_fname,t_test_id,varargin{:});

% Manage VERBLEVEL
if isempty(p.Results.verblevel),
    VERBLEVEL=2; %not global, just local
else
    VERBLEVEL=p.Results.verblevel;
end

%Load GND or GRP struct
if ischar(GND_GRP_or_fname),
    fprintf('Loading GND struct from file %s.\n',GND_GRP_or_fname);
    load(GND_GRP_or_fname,'-MAT');
    if ~exist('GND','var') && ~exist('GRP','var')
        error('File %s does not contain a GND or GRP variable.',GND_GRP_or_fname);
    end
    if exist('GRP','var'),
        GND=GRP; %for simplicity GRP variables are re-named GND
        clear GRP;
    end
    VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL);
else
    GND=GND_GRP_or_fname; %for simplicity GRP and GND variables are re-named GND
    clear GND_GRP_or_fname;
end


if ~strcmpi(GND.t_tests(t_test_id).mean_wind,'yes'),
   error('sig_topo.m only works for significance tests applied to data averaged across multiple time points.'); 
end

n_topos=length(GND.t_tests(t_test_id).used_tpt_ids);
wdth=floor(sqrt(n_topos));
ht=ceil(n_topos/wdth);

bin=GND.t_tests(t_test_id).bin;
chans=GND.t_tests(t_test_id).used_chan_ids;
if isempty(p.Results.fig_id)
    fig_h=figure;
else
    fig_h=figure(p.Results.fig_id); clf;
end
set(fig_h,'name',['Bin ' int2str(bin) ' [' GND.bin_info(bin).bindesc ']'],'paperpositionmode','auto');
%setting paperpositionmode to 'auto' means that if the figure is manually
%resized, the printed version of the figure will reflect whatever the
%shown size was (at the time of printing)

if strcmpi(p.Results.units,'uV'),
    units='\muV';
    %Compute mean voltage in each time window
    mns=zeros(length(chans),n_topos);
    for w=1:n_topos,
        mns(:,w)=mean(GND.grands(chans,GND.t_tests(t_test_id).used_tpt_ids{w},bin),2);
    end
    if isempty(p.Results.one_scale) || (p.Results.one_scale<=0),
        one_scale=0;
        maplimits='absmax';
    else
        one_scale=1;
        mx=max(max(abs(mns)));
        maplimits=[-1 1]*mx;
    end
else
    units='t';
    if isempty(p.Results.one_scale) || (p.Results.one_scale>0),
        one_scale=1;
        mx=max(max(abs(GND.t_tests(t_test_id).data_t)));
        maplimits=[-1 1]*mx;
    else
        one_scale=0;
        maplimits='absmax';
    end
end

if VERBLEVEL,
    if isnan(GND.t_tests(t_test_id).n_perm),
        fprintf('q level of critical t-scores: %f., FDR method: %s\n', ...
            GND.t_tests(t_test_id).desired_alphaORq,GND.t_tests(t_test_id).mult_comp_method);
    else
        fprintf('Estimated alpha level of critical t-scores: %f.\n', ...
            GND.t_tests(t_test_id).estimated_alpha);
    end
end

for a=1:n_topos,
    subplot(ht,wdth,a);
    time_wind=GND.t_tests(t_test_id).time_wind(a,:);
    if strcmpi(p.Results.units,'uV'),
        mn=mns(:,a);
    else
        mn=GND.t_tests(t_test_id).data_t(:,a);
    end
    
    if isnan(GND.t_tests(t_test_id).n_perm),
        sig_chans=find(GND.t_tests(t_test_id).fdr_rej(:,a));
    else
        sig_chans=find(GND.t_tests(t_test_id).adj_pval(:,a)<GND.t_tests(t_test_id).estimated_alpha);
    end
    topoplotMK(mn,GND.chanlocs(chans),'emarker2',{sig_chans,'o',[1 1 1],4}, ...
        'maplimits',maplimits);
    hh=title(sprintf('%d-%d ms',time_wind(1),time_wind(2)));
    if ~one_scale,
        hcb=colorbar;
        hy=ylabel(hcb,units);
        set(hy,'rotation',0,'verticalalignment','middle','position',[4.9+n_topos*.2 0.03 1.0001]);
    end
end

if one_scale,
    hcb=cbar_nudge('vert',0,maplimits);
    hy=ylabel(hcb,units);
    set(hy,'rotation',0,'verticalalignment','middle');
end

% Plot title of entire figure
if p.Results.title_on,
    h=textsc(['Bin ' int2str(bin) ': ' GND.bin_info(bin).bindesc],'title');
    set(h,'fontsize',14,'fontweight','bold');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% cbar_nudge (plots colorbar slightly to left of cbar.m
% search for "DG change" to find the line I edited
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handle]=cbar_nudge(arg,colors,minmax, grad)

if nargin < 2
  colors = 0;
end
posscale = 'off';
if nargin < 1
  arg = 'vert';
  ax = [];
else
  if isempty(arg)
    arg = 0;
  end
  if arg(1) == 0
    ax = [];
    arg = 'vert';
  elseif strcmpi(arg, 'pos')
    ax = [];
    arg = 'vert';
    posscale = 'on';
  else      
    if isstr(arg)
      ax = [];
    else
      ax = arg;
      arg = [];
    end
  end
end

if nargin>2
  if size(minmax,1) ~= 1 | size(minmax,2) ~= 2
    help cbar
    fprintf('cbar() : minmax arg must be [min,max]\n');
    return
  end
end
if nargin < 4
    grad = 5;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Choose colorbar position
%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (length(colors) == 1) & (colors == 0)
  t = caxis;
else
  t = [0 1];
end
if ~isempty(arg)
  if strcmp(arg,'vert')  
    cax = gca;
    pos = get(cax,'Position');
    stripe = 0.04; 
    edge = 0.01;
    space = -.02;

    set(cax,'Position',[pos(1) pos(2) pos(3) pos(4)])
    rect = [pos(1)+pos(3)+space pos(2) stripe*pos(3) pos(4)];
    ax = axes('Position', rect);
  elseif strcmp(arg,'horiz')
    cax = gca;
    pos = get(cax,'Position');
    stripe = 0.075; 
    space = .1;  
    set(cax,'Position',...
        [pos(1) pos(2)+(stripe+space)*pos(4) pos(3) (1-stripe-space)*pos(4)])
    rect = [pos(1) pos(2) pos(3) stripe*pos(4)];
    ax = axes('Position', rect);
  end
else
  pos = get(ax,'Position');
  if pos(3) > pos(4)
    arg = 'horiz';
  else
    arg = 'vert';
  end
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw colorbar using image()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

map = colormap;
n = size(map,1);

if length(colors) == 1
  if strcmp(arg,'vert')
      if strcmpi(posscale, 'on')
          image([0 1],[0 t(2)],[ceil(n/2):n-colors]');
      else
          image([0 1],t,[1:n-colors]');
      end;
      set(ax,'xticklabelmode','manual')
      set(ax,'xticklabel',[],'YAxisLocation','right')
      
  else
    image(t,[0 1],[1:n-colors]);
    set(ax,'yticklabelmode','manual')
    set(ax,'yticklabel',[],'YAxisLocation','right')
  end
  set(ax,'Ydir','normal','YAxisLocation','right')

else % length > 1

  if max(colors) > n
    error('Color vector excedes size of colormap')
  end
  if strcmp(arg,'vert')
    image([0 1],t,[colors]');
    set(ax,'xticklabelmode','manual')
    set(ax,'xticklabel',[])
  else
    image([0 1],t,[colors]);
    set(ax,'yticklabelmode','manual')
    set(ax,'yticklabel',[],'YAxisLocation','right')
  end  
  set(ax,'Ydir','normal','YAxisLocation','right')
end

%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjust cbar ticklabels
%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin > 2 
  if strcmp(arg,'vert')
      Cax = get(ax,'Ylim');
  else
      Cax = get(ax,'Xlim');
  end;
  CBTicks = [Cax(1):(Cax(2)-Cax(1))/(grad-1):Cax(2)]; % caxis tick positions
  CBLabels = [minmax(1):(minmax(2)-minmax(1))/(grad-1):minmax(2)]; % tick labels
  
  dec = floor(log10(max(abs(minmax)))); % decade of largest abs value
  CBLabels = ([minmax]* [ linspace(1,0, grad);linspace(0, 1, grad)]);
  if dec<1
    CBLabels = round(CBLabels*10^(1-dec))*10^(dec-1);
  elseif dec == 1
    CBLabels = round(CBLabels*10^(2-dec))*10^(dec-2);
  else
    CBLabels = round(CBLabels);
  end

  if strcmp(arg,'vert')
    set(ax,'Ytick',CBTicks);
    set(ax,'Yticklabel',CBLabels);
  else
    set(ax,'Xtick',CBTicks);
    set(ax,'Xticklabel',CBLabels);
  end
end
handle = ax;

%%%%%%%%%%%%%%%%%%
% Adjust cbar tag
%%%%%%%%%%%%%%%%%%

set(ax,'tag','cbar')