function [dtf, dtfvar, n] = ft_connectivity_dtf(inputdata, varargin)

% FT_CONNECTIVITY_DTF computes the directed transfer function.
%
% Use as
%   [d, v, n] = ft_connectivity_dtf(inputdata, ...)
%
% The input should be a spectral transfer matrix organized as
%   Nrpt x Nchan x Nchan x Nfreq (x Ntime)
% where Nrpt can be 1.
%
% The output represents
%   d = partial directed coherence matrix Nchan x Nchan x Nfreq (x Ntime).
%       If multiple observations in the input, the average is returned.
%   v = variance of d across observations.
%   n = number of observations.
%
% Typically, nrpt should be 1 where the spectral transfer matrix is computed across
% observations. When nrpt>1 and hasjack=true, the input is assumed to contain the
% leave-one-out estimates of the spectral transfer matrix, thus a more reliable
% estimate of the relevant quantities.
%
% Additional optional input arguments come as key-value pairs:
%   'hasjack'   = boolean, specifying whether the input contains leave-one-outs,
%                 required for correct variance estimate (default = false)
%   'crsspctrm' = matrix containing the cross-spectral density. If this
%                 matrix is defined, the function returns the ddtf, which
%                 requires an estimation of partial coherence from this matrix.
%   'invfun'    = 'inv' (default) or 'pinv', the function used to invert the
%                 crsspctrm matrix to obtain the partial coherence. Pinv is
%                 useful if the data are poorly-conditioned.
%   'feedback'  = 'none', 'text', 'textbar', 'dial', 'etf', 'gui' type of feedback showing progress of computation, see FT_PROGRESS
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS

% Copyright (C) 2009-2017, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

hasjack   = ft_getopt(varargin, 'hasjack', 0);
powindx   = ft_getopt(varargin, 'powindx');
feedback  = ft_getopt(varargin, 'feedback', 'none');
crsspctrm = ft_getopt(varargin, 'crsspctrm');
invfun    = ft_getopt(varargin, 'invfun', 'inv');

switch invfun
  case {'inv' 'pinv'}
    invfun = str2func(invfun);
  otherwise
    ft_error('unknown specification of inverse function for the transfer matrix');
end

if ~isempty(powindx)
  % this error message is rather uninformative, but is kept for now for
  % backward compatibility reasons (i.e. it might exist when called from
  % ft_connectivityanalysis
  ft_error('linearly indexed data for DTF computation is at the moment not supported');
end

siz    = [size(inputdata) 1];
n      = siz(1);
outsum = zeros(siz(2:end));
outssq = zeros(siz(2:end));

% check the crsspctrm; if it is present, compute the partial coherence
if ~isempty(crsspctrm)
  assert(isequal(size(crsspctrm),size(inputdata)), 'the input data should be of the same size as the crsspctrm');
  fprintf('computing dDTF in the presence of a crsspctrm\n');
  
  % the crsspctrm allows for the partial coherence to be computed
  pdim   = prod(siz(4:end));
  tmpcrs = reshape(crsspctrm, [siz(1:3) pdim]);
  ft_progress('init', feedback, 'computing partial coherence...');
  for k = 1:n
    ft_progress(k/n, 'computing partial coherence for replicate %d from %d\n', k, n);
    tmp = reshape(tmpcrs(k,:,:,:), [siz(2:3) pdim]);
    for m = 1:pdim
      tmp(:,:,m) = invfun(tmp(:,:,m));
      tmp(:,:,m) = abs(tmp(:,:,m))./sqrt(abs(diag(tmp(:,:,m))*diag(tmp(:,:,m))'));
    end
    tmpcrs(k,:,:,:) = tmp;
  end
  ft_progress('close');
  crsspctrm = reshape(tmpcrs, siz);
end

% data should be represented as chan_chan_therest
for j = 1:n
  tmph   = reshape(inputdata(j,:,:,:,:), siz(2:end));
  if isempty(crsspctrm)
    % plain DTF
    den    = sum(abs(tmph).^2,2);
    tmpdtf = abs(tmph)./sqrt(repmat(den, [1 siz(2) 1 1 1]));
  else
    % dDTF
    tmpc   = reshape(crsspctrm(j,:,:,:,:), siz(2:end));
    den    = sum(sum(abs(tmph).^2,3),2);
    tmpdtf = abs(tmph)./sqrt(repmat(den, [1 siz(2) siz(4) 1 1 1]));
    tmpdtf = tmpdtf.*tmpc;
  end
  outsum = outsum + tmpdtf;
  outssq = outssq + tmpdtf.^2;
end
dtf = outsum./n;

if n>1 % FIXME this is strictly only true for jackknife, otherwise other bias is needed
  if hasjack
    bias = (n - 1).^2;
  else
    bias = 1;
  end
  dtfvar = bias.*(outssq - (outsum.^2)/n)./(n-1);
else
  dtfvar = [];
end

% this is to achieve the same convention for all connectivity metrics:
% row -> column
for k = 1:prod(siz(4:end))
  dtf(:,:,k)    = transpose(dtf(:,:,k));
  if ~isempty(dtfvar)
    dtfvar(:,:,k) = transpose(dtfvar(:,:,k));
  end
end
