function [freq] = ft_spiketriggeredspectrum_stat(cfg,spike)

% FT_SPIKETRIGGEREDSPECTRUM_STAT computes phase-locking statistics for spike-LFP
% phases. These contain the PPC statistics according to Vinck et al. 2010 (Neuroimage)
% and Vinck et al. 2011 (Journal of Computational Neuroscience).
%
% Use as:
%   [stat] = ft_spiketriggeredspectrum_stat(cfg, spike)
%
% Inputs:
%   SPIKE should be a structure as obtained from from the FT_SPIKETRIGGEREDSPECTRUM function.
%
% Configurations (cfg) 
%
% cfg.method  = string, indicating which statistic to compute. Can be:
%     -'plv' : phase-locking value, computes the resultant length over spike
%              phases. More spikes -> lower value (bias).
%     -'ang' : computes the angular mean of the spike phases.
%     -'ral' : computes the rayleigh p-value.
%     -'ppc0': computes the pairwise-phase consistency across all available
%              spike pairs (Vinck et al., 2010).
%     -'ppc1': computes the pairwise-phase consistency across all available
%              spike pairs with exclusion of spike pairs in the same trial.
%              This avoids history effects within spike-trains to influence
%              phase lock statistics.
%     -'ppc2': computes the PPC across all spike pairs with exclusion of
%              spike pairs in the same trial, but applies a normalization
%              for every set of trials. This estimator has more variance but
%              is more robust against dependencies between spike phase and
%              spike count.
%         
% cfg.timwin  = double or 'all' (default)
%   - double: indicates we compute statistic with a
%            sliding window of cfg.timwin, i.e. time-resolved analysis.
%   - 'all': we compute statistic over all time-points,
%            i.e. in non-time resolved fashion.
%
% cfg.winstepsize  = double, stepsize of sliding window in seconds. For
%   example if cfg.winstepsize = 0.1, we compute stat every other 100 ms.
%
% cfg.channel      = Nx1 cell-array or numerical array with selection of
%   channels (default = 'all'),See CHANNELSELECTION for details
%
% cfg.spikechannel = label of ONE unit, according to FT_CHANNELSELECTION
%
% cfg.spikesel     = 'all' (default) or numerical or logical selection of spikes.
%
% cfg.foi          = 'all' or numerical vector that specifies a subset of
%   frequencies in Hz, e.g. cfg.foi = spike.freq(1:10);                                    
%
% cfg.latency      = [beg end] in sec, or 'maxperiod',  'poststim' or
%  'prestim'.  This determines the start and end of analysis window.
%
% cfg.avgoverchan  = 'weighted', 'unweighted' (default) or 'no'.
%                  This regulates averaging of fourierspectra prior to
%                  computing the statistic.
%  - 'weighted'  : we average across channels by weighting by the LFP power.
%                  This is identical to adding the raw LFP signals in time 
%                  and then taking their FFT.
%  - 'unweighted': we average across channels after normalizing for LFP power. 
%                  This is identical to normalizing LFP signals for 
%                  their power, averaging them, and then taking their FFT.
%  - 'no'        : no weighting is performed, statistic is computed for
%                  every LFP channel.
% cfg.trials       = vector of indices (e.g., 1:2:10),
%                    logical selection of trials (e.g., [1010101010]), or
%                   'all' (default)
%
% Main outputs:
%   stat.nspikes                    =  nChancmb-by-nFreqs-nTimepoints number
%                                      of spikes used to compute stat
%   stat.dimord                     = 'chan_freq_time'
%   stat.labelcmb                   =  nChancmbs cell array with spike vs
%                                      LFP labels
%   stat.(cfg.method)               =  nChancmb-by-nFreqs-nTimepoints  statistic
%   stat.freq                       =  1xnFreqs array of frequencies
%   stat.nspikes                    =  number of spikes used to compute
%
% The output stat structure can be plotted using ft_singleplotTFR or ft_multiplotTFR.

% Copyright (C) 2012, 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_spiketriggeredspectrum_stat.m 10448 2015-06-10 18:01:56Z roboos $

revision = '$Id: ft_spiketriggeredspectrum_stat.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 if the data is of spike format, and convert from old format if required
spike = ft_checkdata(spike,'datatype', 'spike', 'feedback', 'yes');

% get the options
cfg.method         = ft_getopt(cfg, 'method', 'ppc1');
cfg.channel        = ft_getopt(cfg, 'channel', 'all');
cfg.spikechannel   = ft_getopt(cfg, 'spikechannel', spike.label{1});
cfg.latency        = ft_getopt(cfg, 'latency', 'maxperiod');
cfg.spikesel       = ft_getopt(cfg, 'spikesel', 'all');
cfg.avgoverchan    = ft_getopt(cfg, 'avgoverchan', 'no');
cfg.foi            = ft_getopt(cfg, 'foi', 'all');
cfg.trials         = ft_getopt(cfg, 'trials', 'all');
cfg.timwin         = ft_getopt(cfg, 'timwin', 'all');
cfg.winstepsize    = ft_getopt(cfg, 'winstepsize', 0.01);

% ensure that the options are valid
cfg = ft_checkopt(cfg, 'method', 'char', {'ppc0', 'ppc1', 'ppc2', 'ang', 'ral', 'plv'});
cfg = ft_checkopt(cfg, 'foi',{'char', 'double'});
cfg = ft_checkopt(cfg, 'spikechannel',{'cell', 'char', 'double'});
cfg = ft_checkopt(cfg, 'channel', {'cell', 'char', 'double'});
cfg = ft_checkopt(cfg, 'spikesel', {'char', 'logical', 'double'});
cfg = ft_checkopt(cfg, 'avgoverchan', 'char', {'weighted', 'unweighted', 'no'});
cfg = ft_checkopt(cfg, 'latency', {'char', 'doublevector'});
cfg = ft_checkopt(cfg, 'trials', {'char', 'double', 'logical'}); 
cfg = ft_checkopt(cfg, 'timwin', {'double', 'char'}); 
cfg = ft_checkopt(cfg, 'winstepsize', {'double'}); 

cfg = ft_checkconfig(cfg, 'allowed', {'method', 'channel', 'spikechannel', 'latency', 'spikesel', 'avgoverchan', 'foi', 'trials', 'timwin', 'winstepsize'});

% collect channel information
cfg.channel        = ft_channelselection(cfg.channel, spike.lfplabel);
chansel            = match_str(spike.lfplabel, cfg.channel); 

% get the spikechannels
spikelabel       = spike.label;
cfg.spikechannel = ft_channelselection(cfg.spikechannel, spikelabel);
unitsel          = match_str(spikelabel, cfg.spikechannel);
nspikesel        = length(unitsel); % number of spike channels
if nspikesel>1, error('only one unit should be selected for now'); end
if nspikesel==0, error('no unit was selected'); end

% collect frequency information
if strcmp(cfg.foi, 'all'),  
  cfg.foi           = spike.freq(:)'; 
  freqindx          = 1:length(spike.freq);
else
  freqindx  = zeros(1,length(foi));
  for iFreq = 1:length(cfg.foi)
    freqindx(iFreq)         = nearest(spike.freq,cfg.foi(iFreq)); 
  end
end
if length(freqindx)~=length(unique(freqindx)) 
  error('Please select every frequency only once, are you sure you selected in Hz?')
end
nFreqs     = length(freqindx);
cfg.foi    = spike.freq(freqindx); % update the information again

% create the spike selection 
nSpikes = length(spike.trial{unitsel});
if strcmp(cfg.spikesel,'all'), 
  cfg.spikesel = 1:length(spike.trial{unitsel});
elseif islogical(cfg.spikesel)
  cfg.spikesel = find(cfg.spikesel);
end
if ~isempty(cfg.spikesel)
  if max(cfg.spikesel)>nSpikes || length(cfg.spikesel)>nSpikes 
    error('cfg.spikesel must not exceed number of spikes and select every spike just once')
  end
end
  
% select on basis of latency
if strcmp(cfg.latency, 'maxperiod'),
   cfg.latency = [min(spike.trialtime(:)) max(spike.trialtime(:))];
elseif strcmp(cfg.latency,'prestim')
   cfg.latency = [min(spike.trialtime(:)) 0];
elseif strcmp(cfg.latency,'poststim')
   cfg.latency = [0 max(spike.trialtime(:))];
elseif ~isrealvec(cfg.latency) && length(cfg.latency)~=2 
   error('cfg.latency should be "maxperiod", "prestim", "poststim" or 1-by-2 numerical vector');
end
if cfg.latency(1)>=cfg.latency(2), 
   error('cfg.latency should be a vector in ascending order, i.e., cfg.latency(2)>cfg.latency(1)');
end
inWindow = find(spike.time{unitsel}>=cfg.latency(1) & cfg.latency(2)>=spike.time{unitsel});

% selection of the trials
cfg        = trialselection(cfg,spike);

% do the final selection, and select on spike structure
isintrial    = ismember(spike.trial{unitsel}, cfg.trials);
spikesel     = intersect(cfg.spikesel(:),inWindow(:));
spikesel     = intersect(spikesel,find(isintrial));
cfg.spikesel = spikesel; 
spikenum     = length(cfg.spikesel); % number of spikes that were finally selected
if isempty(spikenum), warning('No spikes were selected after applying cfg.latency, cfg.spikesel and cfg.trials'); end
spike.fourierspctrm = spike.fourierspctrm{unitsel}(cfg.spikesel,chansel,freqindx);
spike.time          = spike.time{unitsel}(cfg.spikesel);
spike.trial         = spike.trial{unitsel}(cfg.spikesel);

% average the lfp channels (weighted, unweighted, or not)
if strcmp(cfg.avgoverchan,'unweighted')
  spike.fourierspctrm = spike.fourierspctrm ./ abs(spike.fourierspctrm); % normalize the angles before averaging   
  spike.fourierspctrm = nanmean(spike.fourierspctrm,2); % now rep x 1 x freq
  nChans = 1;
  outlabels = {'avgLFP'}; 
elseif strcmp(cfg.avgoverchan,'no')
  nChans = length(chansel);
  outlabels = cfg.channel;
elseif strcmp(cfg.avgoverchan,'weighted')
  spike.fourierspctrm = nanmean(spike.fourierspctrm,2); % now rep x 1 x freq
  outlabels = {'avgLFP'}; 
  nChans = 1;
end

% normalize the spectrum first
spike.fourierspctrm = spike.fourierspctrm ./ abs(spike.fourierspctrm); % normalize the angles before averaging   
ft_progress('init', 'text',     'Please wait...');
if strcmp(cfg.timwin,'all')
  freq.time = 'all';
  switch cfg.method
    case 'ang'
     out  = angularmean(spike.fourierspctrm);
    case 'plv'
     out  = resultantlength(spike.fourierspctrm);
    case 'ral'
     out  = rayleightest(spike.fourierspctrm); % the rayleigh test
    case 'ppc0'
     out  = ppc(spike.fourierspctrm);
    case {'ppc1', 'ppc2'}
      
      % check the final set of trials present in the spikes
      trials = unique(spike.trial);

      % loop init for PPC 2.0
      [S,SS,dof,dofSS] = deal(zeros(1,nChans,nFreqs));
      nTrials = length(trials);
      if nTrials==1
        warning('computing ppc1 or ppc2 can only be performed with more than 1 trial');
      end
      for iTrial = 1:nTrials % compute the firing rate
          trialNum      = trials(iTrial);
          spikesInTrial = find(spike.trial == trialNum);
          spc           = spike.fourierspctrm(spikesInTrial,:,:);
          ft_progress(iTrial/nTrials, 'Processing trial %d from %d', iTrial, nTrials);              
          % compute PPC 1.0 and 2.0 according to Vinck et al. (2011) using summation per trial
          if strcmp(cfg.method,'ppc2')
            if ~isempty(spc)
              m = nanmean(spc,1); % no problem with NaN
              hasNum = ~isnan(m);
              S(hasNum)  = S(hasNum)  + m(hasNum); % add the imaginary and real components
              SS(hasNum) = SS(hasNum) + m(hasNum).*conj(m(hasNum));
              dof(hasNum) = dof(hasNum) + 1; % dof needs to be kept per frequency
            end
          elseif strcmp(cfg.method,'ppc1')
            if ~isempty(spc)
               n      = sum(~isnan(spc),1);
               m      = nansum(spc,1); 
               hasNum = ~isnan(m);    
               S(hasNum)     = S(hasNum)  + m(hasNum); % add the imaginary and real components
               SS(hasNum)    = SS(hasNum) + m(hasNum).*conj(m(hasNum));
               dof(hasNum)   = dof(hasNum)  + n(hasNum);
               dofSS(hasNum) = dofSS(hasNum) + n(hasNum).^2;
            end                              
          end
      end

      [out]   = deal(NaN(1,nChans,nFreqs));
      hasTrl = dof>1;
      if strcmp(cfg.method,'ppc1')
        out(hasTrl) = (S(hasTrl).*conj(S(hasTrl)) - SS(hasTrl))./(dof(hasTrl).^2 - dofSS(hasTrl));
      elseif strcmp(cfg.method, 'ppc2')
        out(hasTrl) = (S(hasTrl).*conj(S(hasTrl)) - SS(hasTrl))./(dof(hasTrl).*(dof(hasTrl)-1));
      end
  end
  nSpikes = sum(~isnan(spike.fourierspctrm));
else % compute time-resolved spectra of statistic
  
  % make the sampling axis for the window
  bins      = cfg.latency(1):cfg.winstepsize:cfg.latency(2);
  N         = length(bins)-1; % number of bins
  wintime   = 0:cfg.winstepsize:cfg.timwin;
  win       = ones(1,length(wintime));
  freq.time = (bins(2:end)+bins(1:end-1))/2;
  if ~mod(length(win),2), win = [win 1]; end   % make sure the number of samples is uneven.

  out      = NaN(N,nChans,nFreqs);
  nSpikes  = zeros(N,nChans,nFreqs);
  for iChan = 1:nChans
    for iFreq = 1:nFreqs
  
      spctra       = spike.fourierspctrm(:,iChan,iFreq); % values to accumulate
      tm           = spike.time;
      hasnan       = isnan(spctra);    
      tm(hasnan) = [];
      spctra(hasnan) = [];

      % compute the degree of freedom per time bin and the index for every spike
      [dof, indx] = histc(tm, bins); % get the index per spike, and number per bin                          
      if isempty(dof), continue,end    
      toDel       = indx==length(bins) | indx==0; % delete those points that are equal to last output histc or don't fall in
      spctra(toDel) = [];
      indx(toDel)   = [];
      dof(end)    = []; % the last bin is a single point in time, so we delete it

      % compute the sum of spikes per window at every time point
      dof    = dof(:); % force it to be row
      dof    = conv2(dof(:),win(:),'same'); % get the total number of spikes across all trials
      nSpikes(:,iChan,iFreq) = dof;      
      
      switch cfg.method 
        case {'ang', 'plv', 'ppc0', 'ral'}

          % first create a vector with the phases at the samples
          x    = accumarray(indx(:),spctra,[N 1]); % simply the sum of the complex numbers
          
          % then compute the sum of the spectra for every timepoint
          y    = conv2(x(:),win(:),'same');     
          
          % now compute the output statistic
          hasnum  = dof>1;
          hasnum0 = dof>0;          
          if strcmp(cfg.method,'plv')
            out(hasnum,iChan,iFreq) = abs(y(hasnum)./dof(hasnum));
          elseif strcmp(cfg.method,'ral')
            Z = abs(y).^2./dof;
            P = exp(-Z) .* (1 + (2*Z - Z.^2)./(4*dof) -(24.*Z - 123*(Z.^2) + 76*(Z.^3) - 9*(Z.^4))./(288*dof.^2)); %Mardia 1972
            out(hasnum,iChan,iFreq) = P(hasnum);     
          elseif strcmp(cfg.method,'ang')
            out(hasnum0,iChan,iFreq) = angle(y(hasnum0)); %            
          elseif strcmp(cfg.method,'ppc0')
            out(hasnum,iChan,iFreq) = (y(hasnum).*conj(y(hasnum)) - dof(hasnum))./(dof(hasnum).*(dof(hasnum)-1)); % simplest form of ppc
          end

        case {'ppc1', 'ppc2'}
          trials  = unique(spike.trial);
          nTrials = length(trials);
          [S,SS,dofS,dofSS] = deal(zeros(length(bins)-1,1));
          if nTrials==1
             warning('computing ppc1 or ppc2 can only be performed with more than 1 trial');
          end      
          % compute the new ppc versions
          for iTrial = 1:nTrials

            ft_progress(iTrial/nTrials, 'Processing trial %d from %d for freq %d and chan %d', iTrial, nTrials, iFreq, iChan);                
             % select the spectra, time points, and trial numbers again
             trialNum      = trials(iTrial);
             spikesInTrial = find(spike.trial == trialNum);
             if isempty(spikesInTrial), continue,end
             spctraTrial  = spike.fourierspctrm(spikesInTrial,iChan,iFreq);       
             tm           = spike.time;
             hasnan       = isnan(spctraTrial);    
             tm(hasnan) = [];
             spctraTrial(hasnan) = [];

             % bin the spikes and delete spikes out of the selected time
             [dof, indx]     = histc(tm(spikesInTrial), bins); % get the index per spike, and number per bin                          
             dof(end)        = [];
             toDel           = indx==length(bins) | indx==0; % delete those points that are equal to last output histc
             spctraTrial(toDel,:,:) = []; % delete those spikes from fourierspctrm as well
             indx(toDel) = []; % make sure index doesn't contain them
                          
             % first create a vector that sums the spctra
             x    = accumarray(indx(:),spctraTrial,[N 1]); 

             % then compute the moving average sum of this vector
             y    = conv2(x(:),win(:),'same'); % convolution again just means a sum                                  
             d    = conv2(dof(:),win(:),'same'); % get the dof of spikes per trial

             if strcmp(cfg.method,'ppc1')
               S     = S     + y; % add the imaginary and real components
               SS    = SS    + y.*conj(y);
               dofS  = dofS  + d(:);
               dofSS = dofSS + d(:).^2;
             else
               sl             = d(:)>0;
               m              = y./d(:);
               S(sl)          = S(sl)  + m(sl); % add the imaginary and real components
               SS(sl)         = SS(sl) + m(sl).*conj(m(sl));
               dofS(sl)       = dofS(sl) + 1;                     
             end
          end
          
          % we need at least two trials
          hasNum = dofS>1;
          if strcmp(cfg.method,'ppc1')
            out(hasNum,iChan,iFreq) = (S(hasNum).*conj(S(hasNum)) - SS(hasNum))./(dofS(hasNum).^2 - dofSS(hasNum));
          else
            out(hasNum,iChan,iFreq) = (S(hasNum).*conj(S(hasNum)) - SS(hasNum))./(dofS(hasNum).*(dofS(hasNum)-1));
          end
      end
    end
  end
end
ft_progress('close');

% collect the outputs: in labelcmb representation
outparam        = cfg.method;
freq.(outparam) = permute(out,[2 3 1]);
freq.nspikes    = permute(nSpikes,[2 3 1]);    % also cross-unit purposes
freq.labelcmb = cell(nChans,2);
freq.labelcmb(1:nChans,1) = cfg.spikechannel;
for iCmb = 1:nChans
  freq.labelcmb{iCmb,2}   = outlabels{iCmb}; 
end  
freq.freq       = spike.freq(freqindx);
freq.dimord     = 'chancmb_freq_time';

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


function [P] = rayleightest(x)

n = sum(~isnan(x),1);
R = resultantlength(x);    
Z = n.*R.^2;
    
P = exp(-Z).*...
(1 + (2*Z - Z.^2)./(4*n) -(24.*Z - 123*(Z.^2) + 76*(Z.^3) - 9*(Z.^4))./(288*n.^2)); %Mardia 1972

function [resLen] = resultantlength(angles)

n = sum(~isnan(angles),1);
resLen = abs(nansum(angles,1))./n;%calculate the circular variance

function [y] = ppc(crss)

dim = 1;
dof = sum(~isnan(crss),dim);
sinSum = abs(nansum(imag(crss),dim));
cosSum = nansum(real(crss),dim);
y = (cosSum.^2+sinSum.^2 - dof)./(dof.*(dof-1));     
 
function [angMean] = angularmean(angles)

angMean = angle(nanmean(angles,1)); %get the mean angle


function [cfg] = trialselection(cfg,spike)

% get the number of trials or change DATA according to cfg.trials
nTrials = size(spike.trialtime,1);
if  strcmp(cfg.trials, 'all')
  cfg.trials = 1:nTrials;
elseif islogical(cfg.trials)
  cfg.trials = find(cfg.trials);
end
cfg.trials = sort(cfg.trials(:));
if max(cfg.trials)>nTrials, warning('maximum trial number in cfg.trials should not exceed length of DATA.trial')
end
if isempty(cfg.trials), error('No trials were selected');
end

function m = nansum(x,dim)

% Find NaNs and set them to zero
nans = isnan(x);
x(nans) = 0;

if nargin == 1 % let sum deal with figuring out which dimension to use
    % Count up non-NaNs.
    n = sum(~nans);
    n(n==0) = NaN; % prevent divideByZero warnings
    % Sum up non-NaNs, and divide by the number of non-NaNs.
    m = sum(x);
else
    % Count up non-NaNs.
    n = sum(~nans,dim);
    n(n==0) = NaN; % prevent divideByZero warnings
    % Sum up non-NaNs, and divide by the number of non-NaNs.
    m = sum(x,dim);
end






