function [F,q] = lyngby_baf2_test(X, R, arg1, arg2, arg3, arg4, arg5, arg6)

% lyngby_baf2_test      - Barbak Ardekani's F-test with nuisance estimation
%
%       function F = lyngby_baf_test(X, R, 'PropertyName',
%                                          'PropertyValue'...)
%
%       Input:  X   The datamatrix
%               R   The run structure
%               Property:
%                  'Components'     Number of Fourier components used in
%                                   the test
%                  'MaxNuisanceDim' Maximal dimension of the nuisance subspace
%                  'TestLevel'      Test level
%
%       Output: F   F-map
%
%       See also: lyngby_bat_test

% cvs : $Id: lyngby_baf2_test.m,v 1.2 1997/11/24 12:13:16 pto Exp $
%       $Revision: 1.2 $

% Parsing
n = 1;
while n < nargin-1
  eval(sprintf('arg = lower(arg%d);', n)); 
  if strcmp(arg, 'components')
    n = n + 1;
    eval(sprintf('arg = arg%d;', n));
    Components = arg;
    if Components<0 
      error('The number of components should be non-negative')
    end
  elseif strcmp(arg, 'maxnuisancedim')
    n = n + 1;
    eval(sprintf('arg = arg%d;', n));
    MaxNuisanceDim = arg;
    if MaxNuisanceDim<1 
      error('Max dimensions of the nuisance subspace must be positive integer')
    end    
  else
    error('Invalid property');
  end

  n = n + 1;
end

M = Components;


% Find number of scans
[tempv,T] = lyngby_runinfo(R);

if length(T) > 1
  error(['The Ardekani F-tests cannot handle runs with ', ...
        'different lengths.']);
  
end

[NoTime,NoVoxels] = size(X);

msplit = floor(M/2);
omega = 2*pi/T;
in = (0:NoTime-1)';

for k = 1:msplit,
  a(:,k) = cos(k*omega*in);
end

for k = msplit+1:M,
  q = k-msplit;
  a(:,k) = sin(q*omega*in);
end  

lyngby_log(sprintf('Identifying signal subspace'));
Pa = a*inv(a'*a)*a';

I = eye(size(Pa));

F = zeros(1,NoVoxels);

lyngby_log(sprintf('Computing total noise signal'));

reX = Pa*X;
Y = X-reX;

lyngby_log(sprintf('Estimating the signal part (%i/%i)',1,NoVoxels));
for j=1:NoVoxels,
  if lyngby_mod(j,50) == 0,
    lyngby_log(sprintf('Estimating the signal part (%i/%i)',j,NoVoxels));
  end
  Fn(j) = (X(:,j))'*reX(:,j)/Components;
end
clear reX;

C = zeros(NoTime,NoTime);
lyngby_log(sprintf('Computing correlation matrix (%i/%i)',1,NoVoxels));
for ii=1:NoVoxels,
  if lyngby_mod(ii,50) == 0,
    lyngby_log(sprintf('Computing correlation matrix (%i/%i)',ii,NoVoxels));
  end
  C = C + Y(:,ii) * Y(:,ii)';
end
C = C/NoVoxels;

lyngby_log(sprintf('Computing svd of correlation matrix'));
[ub,sb,sv] = svd(C);


lyngby_log(sprintf('Computing AIC (%i/%i)',0,MaxNuisanceDim));
Pe = I-Pa;
remenergy=0;
for k=1:NoVoxels,
  remenergy = remenergy+(X(:,k))'*Pe*X(:,k);
end
AIC_zero = NoVoxels*Components/2*log(remenergy);

for i=1:MaxNuisanceDim,
  lyngby_log(sprintf('Computing AIC (%i/%i)',i,MaxNuisanceDim));
  Pb = ub(:,1:i)*ub(:,1:i)';
  Pe = I-Pa-Pb;
  remenergy(i)=0;
  for k=1:NoVoxels,
    remenergy(i) = remenergy(i)+(X(:,k))'*Pe*X(:,k);
  end
  AIC(i) = NoVoxels*Components/2*log(remenergy(i))+i*(Components+NoVoxels);
end

figure;plot(0:MaxNuisanceDim,[AIC_zero AIC])

title('AIC');
[mv,q]=min(AIC);

if (mv<AIC_zero)
  Pb = ub(:,1:q)*ub(:,1:q)';
else
  Pb = 0;
end
Pe = I-Pa-Pb;

lyngby_log(sprintf('Estimating the F map (%i/%i)',1,NoVoxels));
for j = 1:NoVoxels,
  if lyngby_mod(j,100) == 0,
    lyngby_log(sprintf('Estimating the F map (%i/%i)',j,NoVoxels));
  end
  F(j) = Fn(j) /(X(:,j)'*Pe*X(:,j)) * (NoTime-q-Components);
end







