function [cfg] = ft_spike_plot_psth(cfg, psth)

% FT_SPIKE_PLOT_PSTH makes a bar plot of PSTH structure with error bars.
%
% Use as
%   ft_spike_plot_psth(cfg, psth)
% 
%	Inputs:
%		PSTH typically is a structure from FT_SPIKE_PSTH.
%
% Configurations:
%   cfg.latency          = [begin end] in seconds, 'maxperiod' (default), 'prestim'(t<=0), or
%                          'poststim' (t>=0).
%   cfg.errorbars        = 'no', 'std', 'sem' (default), 'conf95%' (requires statistic toolbox,
%                          according to student-T distribution), 'var'
%   cfg.spikechannel     = string or index of single spike channel to trigger on (default = 1)
%                          Only one spikechannel can be plotted at a time.
%   cfg.ylim             = [min max] or 'auto' (default)
%                          If 'standard', we plot from 0 to 110% of maximum plotted value);
% Outputs:
%	  cfg.hdl.avg              = figure handle for the bar plot, psth average.
%	  cfg.hdl.var              = figure handle for the error lines. 
%
% See also FT_SPIKE_PSTH

% Copyright (C) 2010, Martin Vinck
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ft_spike_plot_psth.m 10448 2015-06-10 18:01:56Z roboos $

revision = '$Id: ft_spike_plot_psth.m 10448 2015-06-10 18:01:56Z roboos $';

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble callinfo
ft_preamble trackconfig

% check the psth datatype
psth = ft_checkdata(psth, 'datatype', 'timelock', 'feedback', 'yes');

% get the default options
cfg.spikechannel = ft_getopt(cfg, 'spikechannel', 'all');
cfg.latency      = ft_getopt(cfg,'latency','maxperiod');
cfg.errorbars    = ft_getopt(cfg,'errorbars', 'sem');
cfg.ylim         = ft_getopt(cfg,'ylim', 'auto');

% ensure that the options are valid
cfg = ft_checkopt(cfg,'spikechannel',{'cell', 'char', 'double'});
cfg = ft_checkopt(cfg,'latency', {'char', 'ascendingdoublebivector'});
cfg = ft_checkopt(cfg,'errorbars', 'char', {'sem', 'std', 'conf95%', 'no', 'var'});
cfg = ft_checkopt(cfg,'ylim', {'char', 'ascendingdoublebivector'});

cfg = ft_checkconfig(cfg, 'allowed', {'spikechannel', 'latency', 'errorbars', 'ylim'});

% select the latencies
if ischar(cfg.latency)
  if strcmp(cfg.latency, 'maxperiod')
    cfg.latency = [min(psth.time) max(psth.time)];
  elseif strcmp(cfg.latency, 'prestim')
    cfg.latency = [min(psth.time) 0];
  elseif strcmp(cfg.latency, 'poststim')
    cfg.latency = [0 max(psth.time)];
  end
end

% check whether the time window fits with the data
if cfg.latency(1) < min(psth.time), cfg.latency(1) = min(psth.time);
  warning('Correcting begin latency of averaging window');
end
if (cfg.latency(2) > max(psth.time)), cfg.latency(2) = max(psth.time);
  warning('Correcting begin latency of averaging window');
end

% get the spikechannels
cfg.spikechannel = ft_channelselection(cfg.spikechannel, psth.label);
spikesel    = match_str(psth.label, cfg.spikechannel);
nUnits      = length(spikesel);
if nUnits~=1, error('You selected more or less than one spikechannel by means of cfg.spikechannel'); end

% select the timepoints within the latencies
timeSel = psth.time>=cfg.latency(1) & psth.time <= cfg.latency(2);
if isempty(timeSel), error('no time points selected, please change cfg.latency or inspect input PSTH');end

% plot the average psth
psthHdl = bar(psth.time(timeSel),psth.avg(spikesel,timeSel),'k');
set(psthHdl,'BarWidth', 1)

% check if fields .dof and .var are correct given cfg.errorbars
if ~strcmp(cfg.errorbars,'no')
  if ~isfield(psth,'var') || ~any(isfinite(psth.var(spikesel,:)))
    error('PSTH should contain field .var with numbers');
  end
  if ~ (strcmp(cfg.errorbars,'std') || strcmp(cfg.errorbars,'var'))
    if ~isfield(psth,'dof') || ~any(isfinite(psth.dof(spikesel,:)))
      error('psth should contain field dof that contains numbers.'); 
    end
  end
end

% add standard error bars to it
x = [psth.time(timeSel);psth.time(timeSel)]; % create x doublets
if strcmp(cfg.errorbars, 'sem')
  err = sqrt(psth.var(spikesel,timeSel)./psth.dof(timeSel));
elseif strcmp(cfg.errorbars, 'std')
  err = sqrt(psth.var(spikesel,timeSel));
elseif strcmp(cfg.errorbars, 'var')
  err = psth.var(spikesel,timeSel);
elseif strcmp(cfg.errorbars, 'conf95%')
  % use a try statement just in case the statistics toolbox doesn't work.
  try
    tCrit = tinv(0.975,psth.dof(timeSel));
    err = tCrit.*sqrt(psth.var(spikesel,timeSel)./psth.dof(timeSel)); % assuming normal distribution, SHOULD BE REPLACED BY STUDENTS-T!
  catch
    err = 0;
  end
else
  err = 0;
end
y = psth.avg(spikesel,timeSel) + err; % create the error values

% plot the errorbars as line plot on top of the bar plot
if ~strcmp(cfg.errorbars,'no')
  hold on
  y = [zeros(1,length(y));y]; % let lines run from 0 to error values
  psthSemHdl = plot(x,y,'k'); % plot the error bars
end

% create labels
xlabel('bin centers (sec)')
try
  ylabel(psth.cfg.outputunit)
catch
  ylabel('firing intensity')
end

% set the axis
if strcmp(cfg.ylim, 'auto'), cfg.ylim = [0 max(y(:))*1.1+eps]; end
if cfg.ylim(2)<max(y(:))
  warning('correcting cfg.ylim(2) as it clips the data'); 
  cfg.ylim(2) = max(y(:))*1.1+eps; 
end

set(gca,'YLim', cfg.ylim, 'XLim', cfg.latency)
set(gca,'TickDir','out', 'Box', 'off')

% store the handles
cfg.hdl.avg    = psthHdl;
if ~strcmp(cfg.errorbars,'no'), cfg.hdl.var = psthSemHdl;  end

% as usual, make sure that panning and zooming does not distort the y limits
set(zoom,'ActionPostCallback',{@mypostcallback,cfg.ylim,cfg.latency});
set(pan,'ActionPostCallback',{@mypostcallback,cfg.ylim,cfg.latency});

% do the general cleanup and bookkeeping at the end of the function
ft_postamble trackconfig
ft_postamble callinfo
ft_postamble previous psth

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = mypostcallback(fig,evd,limY,limX)

% get the x limits and reset them
ax = evd.Axes;
xlim = get(ax(1), 'XLim');
ylim = get(ax(1), 'YLim');
if limX(1)>xlim(1), xlim(1) = limX(1); end
if limX(2)<xlim(2), xlim(2) = limX(2); end
if limY(1)>ylim(1), ylim(1) = limY(1); end
if limY(2)<ylim(2), ylim(2) = limY(2); end

set(ax(1), 'XLim',xlim)
set(ax(1), 'YLim',ylim);

