function [Stat] = ft_spike_rate_orituning(cfg,varargin)

% FT_SPIKE_RATE_ORITUNING computes a model of the firing rate as a function
% of orientation or direction.
%
% Use as
%   [stat] = ft_spike_rate_tuning(cfg, rate1,rate2,...rateN)
%
% Inputs rate should be the output from FT_SPIKE_RATE
% Tip: put singel rate structures in cell array rate and write
% stat = ft_spike_rate_tuning(cfg,rate{:});
%
% Configurations:
%   cfg.stimuli  = should be an 1 x nConditions array of orientations or
%                  directions in radians
%                  varargin{i} corresponds to cfg.stimuli(i)
%   cfg.method   = model to apply, implemented are 'orientation' and 'direction'
%
% Outputs:
%   Stat.ang       = mean angle of orientation / direction (1 x nUnits)
%   Stat.osi       = orientation selectivity index (Womelsdorf et al., 2012,
%                    PNAS), that is resultant length.
%                    if cfg.method = 'orientation', then orientations are
%                    first projected on the unit circle.
%   Stat.di        = direction index, 1 - min/max response

% FIXME: models for contrast etc.

% 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_rate_orituning.m 10447 2015-06-10 17:52:33Z roboos $

revision = '$Id: ft_spike_rate_orituning.m 10447 2015-06-10 17:52:33Z roboos $';

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

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', {'stimuli', 'method'});
cfg = ft_checkopt(cfg,'stimuli','doublevector');
cfg = ft_checkopt(cfg,'method', 'char', {'orientation', 'direction'});

if length(varargin)<2, error('can only compute ori tuning if multiple inputs are specified'); end 
  
% check whether trials were kept in the rate function
for k = 1:length(varargin)
  try
    varargin{k} = rmfield(varargin{k}, 'trial');
  end
end

Tune        = ft_appendtimelock([],varargin{:});
Tune.avg    = Tune.trial;

% get the unique stimuli and the number of directions, these should match
stimuli     = cfg.stimuli(:);
nConditions = length(varargin);
nStimuli    = length(stimuli);
if nConditions~=nStimuli, error('Length of cfg.stimuli should match number of data inputs'); end

nUnits      = size(Tune.avg,2);

% change the directions so it is a circular variable with range 2*pi
if strcmp(cfg.method,'orientation')
  if (max(stimuli)-min(stimuli))>pi
    error('If cfg.tuningtype is "orientation", cfg.stimuli should have range of pi');
  end
  stimuli = stimuli*2; % convert to make it circular (see Womelsdorf et al. 2012, PNAS).
  
  if (max(stimuli)-min(stimuli))<0.5*pi
    warning('Orientations have a range < 1/2 pi. Are you sure this is correct?. Stats will be biased');
  end
  
  % compute the directionality index
  Stat.di = 1 - min(Tune.avg)./max(Tune.avg);
  
  % transform the data into complex numbers to compute resultant length
  z = exp(1i*stimuli(:)*ones(1,nUnits));
  sumAvg = sum(Tune.avg);
  z = Tune.avg.*z./(sumAvg(ones(nStimuli,1),:));
  Stat.osi  = abs(sum(z));
  
  % make preferred angle modulo 2pi, convert it back to range pi and convert to rad or deg
  prefAngle           = angle(nansum(z));
  prefAngle           = mod(prefAngle,2*pi);
  Stat.ang      = prefAngle/2;
elseif strcmp(cfg.method,'direction')
  if (max(stimuli)-min(stimuli))>2*pi
    error('Directions has a range > 2*pi. Are you sure this is radians and not degrees?')
  end
  if (max(stimuli)-min(stimuli))<0.5*pi
    warning('Directions have a range < 1/2 pi. Are you sure this is correct?');
  end
  
  % compute the directionality index
  Stat.di = 1 - min(Tune.avg)./max(Tune.avg);
  
  % transform the data into complex numbers to compute resultant length
  z = exp(1i*stimuli(:)*ones(1,nUnits));
  sumAvg = sum(Tune.avg);
  z = Tune.avg.*z./(sumAvg(ones(nStimuli,1),:));
  Stat.osi  = abs(sum(z));
  
  % make preferred angle modulo 2pi, convert it back to range pi and convert to rad or deg
  prefAngle      = angle(nansum(z));
  Stat.ang = mod(prefAngle,2*pi);
end

Stat.label   = Tune.label;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble trackconfig
ft_postamble callinfo
ft_postamble previous Tune
ft_postamble history Stat

