function [F,df,xX,xCon,beta,V] = spm_reml_ancova(y,P,Fc)
% Classical inference for [hierarchial] linear models using ReML
% FORMAT [F,df,xX,xCon,beta,V] = spm_reml_ancova(y,P,Fc);
%
% y       - (n x 1)     response variable
% P{i}.X  - (n x m)     ith level design matrix i.e:
% P{i}.C  - {q}(n x n)  ith level contraints on the form of Cov{e{i}}
% Fc      - (m x q)     contrast matrix for the last level
%
% F     -  T or F values
% df    -  degrees of freedom
% beta  -  parameter estimates
% xX    -  design matrix structure
% xCon  -  contrast structure
%___________________________________________________________________________
%
% spm_ancova uses a General Linear Model of the form:
%
%                            y = X{1}*b{1} + e{1}
%                         b{1} = X{2}*b{2} + e{2}
%                                 ...
%
%                     b{n - 1} = X{n}*b{n} + e{n}
%
% e{n} ~ N{0,Ce{n}} 
%
% An F ratio is formed using OLS estimators of the parameters and ReML 
% estimators of the hyperparamters.
% If Fc has only one column a T statistic is returned,
%___________________________________________________________________________
% @(#)spm_reml_ancova.m	2.1 Karl Friston 02/01/23


% get ReML hyperparameter estimates
%---------------------------------------------------------------------------
[C P] = spm_PEB(y,P);
L     = P{1}.h;			% ReML Hyperparemter estimates
W     = P{1}.W;			% ReML precisions = ddln(p(y|h)/dhdh

% check covariance constraints - assume i.i.d. errors conforming to X{i}
%---------------------------------------------------------------------------
for i = 1:length(P)
	if ~isfield(P{i},'C')
		[n m] = size(P{i}.X);
		if i == 1
			P{i}.C            = {speye(n,n)};
		else
			for j = 1:m
				k         = find(P{i}.X(:,j));
				P{i}.C{j} = sparse(k,k,1,n,n);
			end
		end
	end
end

% form non-hierachical model
%---------------------------------------------------------------------------
X     = 1;
Q     = {};
for i = 1:length(P)
	
	% covariance components
	%-------------------------------------------------------------------
	for j = 1:length(P{i}.C)
		Q{end + 1} = X*P{i}.C{j}*X';
	end

	% design matrix
	%-------------------------------------------------------------------
	X   = X*P{i}.X;
end



% create design matrix structure and get pseudoinverse
%---------------------------------------------------------------------------
xX    = spm_sp('Set',X);
xX.pX = spm_sp('x-',xX);

% OLS parameter  estimates
%---------------------------------------------------------------------------
beta  = xX.pX*y;

% contrast
%---------------------------------------------------------------------------
xCon  = spm_FcUtil('Set','','F','c',Fc,xX);
h     = spm_FcUtil('Hsqr',xCon,xX);
X1o   = spm_FcUtil('X1o',xCon,xX);


% Note tr{MQ} = tr{h*xX.pX*Q{i}*xX.pX'*h'} because
% M = R0 - R = X1o*pinv(X1o) = xX.pX'*h'*h*xX.pX
%---------------------------------------------------------------------------
V     = sparse(0);
T     = sparse(1,length(Q));
for i = 1:length(Q);
	V     = V + L(i)*Q{i};
	T(i)  = trace(h*xX.pX*Q{i}*xX.pX'*h');
end
V     = V*length(V)/trace(V);
TL    = T*L;


% degrees of freedom
%---------------------------------------------------------------------------
[trMV trMVMV] = spm_SpUtil('trMV',X1o,V);
df            = [trMV^2/trMVMV (TL)^2/(T*inv(W)*T')];

if size(Fc,2) == 1

	% T statistics
	%-------------------------------------------------------------------
	F     = h*beta./sqrt(TL);
else
	% F statistics
	%-------------------------------------------------------------------
	F     = sum((h*beta).^2)./TL;

end 
