function [s, cfg] = ft_statfun_depsamplesFmultivariate(cfg, dat, design)

% FT_STATFUN_DEPSAMPLESFMULTIVARIATE calculates the MANOVA dependent samples
% F-statistic on the biological data in dat (the dependent variable), using the
% information on the independent variable (ivar) in design.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option:
%   cfg.statistic = 'ft_statfun_depsamplesFmultivariate'
%
% You can specify the following configuration options:
%   cfg.contrastcoefs  = matrix of contrast coefficients determining the
%                        effect being tested. The number of columns of this
%                        matrix has to be equal to the number of conditions.
%                        The default is a matrix that specifies the
%                        main effect of the independent variable. This matrix
%                        has size [(ncond-1),ncond].
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or cfg.computeprob='yes':
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1)
%
% The experimental design is specified as:
%   cfg.ivar  = independent variable, row number of the design that contains the labels of the conditions to be compared (default=1)
%   cfg.uvar  = unit variable, row number of design that contains the labels of the units-of-observation, i.e. subjects or trials (default=2)
%
% The labels for the independent variable should be specified as numbers ranging
% from 1 to the number of conditions. The labels for the unit of observation should
% be integers ranging from 1 to the total number of observations (subjects or trials).
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

% Copyright (C) 2006, Eric Maris
%
% 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$

% set the defaults
cfg.computestat    = ft_getopt(cfg, 'computestat', 'yes');
cfg.computecritval = ft_getopt(cfg, 'computecritval', 'no');
cfg.computeprob    = ft_getopt(cfg, 'computeprob', 'no');
cfg.alpha          = ft_getopt(cfg, 'alpha', 0.05);
cfg.tail           = ft_getopt(cfg, 'tail', 1);
cfg.ivar           = ft_getopt(cfg, 'ivar', 1);
cfg.uvar           = ft_getopt(cfg, 'uvar', 2);

nconds = length(unique(design(cfg.ivar,:)));

if ~isfield(cfg,'contrastcoefs')
    % specify the default contrast coefficient matrix
    ncontrasts = nconds-1;
    cfg.contrastcoefs = zeros(ncontrasts,nconds);
    cfg.contrastcoefs(:,1) = 1;
    for contrastindx=1:ncontrasts
        cfg.contrastcoefs(contrastindx,contrastindx+1)=-1;
    end
else
    ncontrasts = size(cfg.contrastcoefs,1);
end

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
    ft_error('P-values can only be calculated if the test statistics are calculated.');
end
if ~isfield(cfg,'uvar') || isempty(cfg.uvar)
    ft_error('uvar must be specified for dependent samples statistics');
end

% perform some checks on the design
nuospercond=zeros(nconds,1);
for condindx=1:nconds
    nuospercond(condindx)=length(find(design(cfg.ivar,:)==condindx));
end
if sum(nuospercond)<size(design,2) || any(nuospercond~=nuospercond(1))
  ft_error('Invalid specification of the design array.');
end
nunits = max(design(cfg.uvar,:));
dfdenom = nunits - ncontrasts;
if dfdenom<1
    ft_error('The data must contain more units-of-observation (usually subjects) than the number of contrasts.')
end
nrepl=nunits*nconds;
if (nrepl~=sum(nuospercond)) || (nrepl~=size(dat,2))
  ft_error('Invalid specification of the design array.');
end
nsmpls = size(dat,1);

if strcmp(cfg.computestat,'yes')
% compute the statistic
    % store the positions of the condition labels nunits-by-nconds array
    poslabelsperunit=zeros(nunits,nconds);
    for condindx=1:nconds
        poslabel=find(design(cfg.ivar,:)==condindx);
        [dum,i]=sort(design(cfg.uvar,poslabel),'ascend');
        poslabelsperunit(:,condindx)=poslabel(i);
    end
    % reshape poslabelsperunit into a row vector that contains the
    % replications of the first condition on the first nunits positions,
    % the replications of the second condition on the second nunits
    % positions, etc.
    poslabelsperunit=reshape(poslabelsperunit,1,nrepl);
    s.stat=zeros(nsmpls,1);
    for smplindx=1:nsmpls
        datonesmpl=reshape(dat(smplindx,poslabelsperunit),nunits,nconds);
        contrasts=datonesmpl*cfg.contrastcoefs';
        contrastavg=mean(contrasts,1);
        dev=contrasts-repmat(contrastavg,nunits,1);
        covmat=(dev'*dev)/(nunits-1);
        s.stat(smplindx)=nunits*contrastavg*inv(covmat)*contrastavg';
    end
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.dfnum = ncontrasts;
  s.dfdenom = nunits - ncontrasts;
  if cfg.tail==-1
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a left tail critical value.');
  end
  if cfg.tail==0
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a two-sided critical value.');
  end
  if cfg.tail==1
    s.critval = ((nunits-1).*ncontrasts./(nunits-ncontrasts)).*finv(1-cfg.alpha,s.dfnum,s.dfdenom);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.dfnum = ncontrasts;
  s.dfdenom = nunits - ncontrasts;
  if cfg.tail==-1
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a left tail p-value.');
  end
  if cfg.tail==0
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a two-sided p-value.');
  end
  if cfg.tail==1
    scaledstat = ((nunits-ncontrasts)./((nunits-1).*ncontrasts)).*s.stat;
    s.prob = 1-fcdf(scaledstat,s.dfnum,s.dfdenom);
  end
end
