function  Result = lyngby_lzit_main(Pa, X, arg1, arg2, arg3, arg4, ...
    arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, ...
    arg15);

% lyngby_lzit_main     - Main function for Lange-Zeger estimation
%
%       function Result = lyngby_lzit_main(P, Y, 'PropertyName',
%          'PropertyValue'); 
%
%       Input:    P   Paradigm
%                 Y   Normalized data matrix
%
%       Property: Iterations   { 90 } The number iterations
%                 MinChange    { 1e-4 } Stopcriterion for the
%                              iteration
%                 StepSize     { 1 } Gradient descent stepsize
%                 Theta1Init   { 10 } Initialization of theta1.
%                 Theta2Init   { 2 } Initialization of theta2.
%
%       Output:   Result    Estimated parameters in matrix 3xN 
%                           [beta; theta1; theta2]'
%
%       Fitting of a gamma kernel for each voxel in the Fourier domain.
%       This function makes some initialization before iterating: 
%       for every voxel calling lyngby_lzit_algo to do the actual fitting.
%
%       Ref: Lange, Zeger (1997), Roy Stat Soc, Appl Stat. 44.
%
%       See also LYNGBY, LYNGBY_LZGS_MAIN, LYNGBY_LZIT_ALGO,
%                LYNGBY_LZITPICKF, LYNGBY_LZIT_FTLAMB,
%                LYNGBY_LZIT_THOPT, LYNGBY_LZIT_NOISE,
%                LYNGBY_LZIT_BETA. 
%
% $Id: lyngby_lzit_main.m,v 1.8 2003/02/12 20:06:49 fnielsen Exp $



    lyngby_log('Initializing Lange-Zeger');

 
    % Default values
    iterations = 90;
    minChange = 1e-4;
    stepSize = 1;
    theta1Init = 10;
    theta2Init = 2;
    
    flag = 1;  % [ {1} | 2 ]  !!!
    % 1: first component along (DC-component)
    % 2: first component not along.
    

    % The number of cycles in the paradigm determined by a look on the
    % fft of the paradigm
    fftP = fft(Pa);
    cycles = min(find(abs(fftP).*(abs(fftP)>eps*10^3))) - 1;
    if cycles < 1
      error('The paradigm should be centered.');
    end 


    % Parse properties
    n = 1;
    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 
      if strcmp(arg, 'theta1init')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	theta1Init = arg;
	if theta1Init<0 
	  error('Theta 1 should be a larger than 0.')
	end
      elseif strcmp(arg, 'theta2init')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	theta2Init = arg;
	if theta2Init<0 
	  error('Theta 2 should be a larger than 0.')
	end
      elseif strcmp(arg, 'iterations')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	iterations = arg;
	if iterations<1 
	  error('The number of iterations should be a value larger than 1.')
	end
      elseif strcmp(arg, 'minchange')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	minChange = arg;
	if minChange <= 0
	  error('The minimum change should be a value larger than 0.')
	end
      elseif strcmp(arg, 'stepsize')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	stepSize = arg;
	if stepSize<0
	  error('The step size should be a value larger than 1.')
	end
      else
	error(sprintf('Invalid property %s (lyngby_lzit_main)',arg));
      end
      n = n + 1;
    end
    

    [numTime, numVoxels] = size(X);

    % Extract non-zero frequency components 
    dPi = lyngby_lzit_pickf(fftP,cycles);

    E=[];
    lyngby_log(sprintf('Computing iterative LZ-map: %d / %d', 1, size(X,2)));

    for v = 1:numVoxels
      if lyngby_mod(v,5)==0,
	lyngby_log(sprintf('Computing iterative LZ-map: %d / %d', v, size(X,2)));
      end;

      x = X(1:size(X,1), v);
      
      dX = fft(x(1:numTime));
      
      fhats = 0; 
      dXi = lyngby_lzit_pickf(dX, cycles);
      if flag == 2
	dYi(1,1) = [0];
      end;
      [beta, t1,t2,d,Q_all] = lyngby_lzit_algo(...
	  numTime, ...
	  cycles, ...
	  stepSize, ...
	  minChange, ...
	  iterations, ...
	  fhats, ...
	  dPi, ...
	  dXi, ...
	  theta1Init,...
	  theta2Init);
      
      % Cheap trick to test whether convergence was met....
      if (abs(beta) > 10 * std(X(:,v)))
	E = [E [0 1 1 0]'];  
      else
	E = [E [beta t1 t2 d]'];  
      end
    end; 

    lyngby_log('Iterative LZ-map completed!');

    Result = E(1:3,:);









