function [F, aic, Ub] = lyngby_baf2_main(X, R, arg1, arg2, arg3, arg4, ...
    arg5, arg6)


% lyngby_baf2_main     - Ardekani's F-test with nuisance estimation
%
%       function F = lyngby_baf2_main(X, R, ...
%          'PropertyName', 'PropertyValue')
%
%       Input:  X    The datamatrix
%               R    The run structure
%               Property:
%                 Components      {8} Number of Fourier components
%                                 used in the test
%                 MaxNuisanceDim  {10} Maximal dimension of the
%                                 nuisance subspace 
%
%       Output: F    F-map
%               aic  AIC for 0..MaxNuisanceDim models
%               Ub   Nuisance sequences
%
%       The function performs an F-test on periodic data with the
%       period define by the run structure 'R'. A global subspace
%       orthogonal to the signal subspace is extracted from the
%       residual.
%
%       Ref: Ardekani BA, Kershaw J, Kashikura K, Kanno I. Activation
%            detection in functional MRI using subspace modeling and 
%            maximum likelihood estimation. IEEE Trans Med Imaging
%            1999 Feb, 18(2):101-114 
% 
%       See also LYNGBY, LYNGBY_BAF_TEST, LYNGBY_BAT_TEST.
%
% $Id: lyngby_baf2_main.m,v 1.5 2001/07/31 15:10:04 fnielsen Exp $


    % Sizes
    [NoTime, NoVoxels] = size(X);
    [rX, cX] = size(X);
    [rR, cR] = size(R);
    
    % Check arguments
    if rX ~= rR 
      error('X and R should have the same number of rows');
    end
    
    % Default values
    comp = 8;
    maxNuisanceDim = 10;
    testLevel = 0.001;
    
    % 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));
	comp = arg;
	if comp < 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

    % 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

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

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

    for k = msplit+1:comp,
      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'));
    if 0
      % This is rather slow
      reX = Pa * X;
    else
      % Faster
      reX = a * (inv(a'*a) * (a'* X));
    end
      
    Y = X - reX;

    
    lyngby_log(sprintf('Estimating the signal part: %i / %i',1,NoVoxels));
    Fn = zeros(NoVoxels,1);     % Preallocation of Fn
    for j = 1:NoVoxels,
      if lyngby_mod(j, 100) == 0,
	lyngby_log(sprintf('Estimating the signal part: %i / %i',j,NoVoxels));
      end
      Fn(j) = (X(:,j))' * reX(:,j) / comp;
    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,100) == 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*comp/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*comp/2*log(remenergy(i))+i*(comp+NoVoxels);
    end

    aic = [ AIC_zero AIC ];
    Ub = Ub(:,1:maxNuisanceDim);
    
    [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-comp);
    end










