function [E, Y, U] = lyngby_nnr_main(x, T, arg1, arg2, arg3, arg4, ...
    arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, ...
    arg15, arg16, arg17, arg18, arg19, arg20)

% lyngby_nnr_main      - Main function for Neural network Regression
%
%	function [E, Y, U] = lyngby_nnr_main(x, T, PropertyName, ...
%	                         PropertyValue)
%
%       Input:  x   Neural network input (the paradigm)
%               T   Target output (the datamatrix)
%               PropertyName:
%                  'GenOptim'      { {Free} | EarlyStop | 
%                                  HiddenUnitsEarlyStop | Pruning |
%                                  Pruning1DRegGridSearch |
%                                  Pruning2DRegGridSearch ]
%                                  Generalization optimization
%                  'Validation'    [ {SingleBlocked}]
%                                  How to compute the generalization
%                  'Lag'           { 7 } Time lag.
%                  'HiddenUnits'   { 3 } Number of hidden units, not
%                                  counting the threshold unit
%                  'Reg'           { 0.001 } Regularization parameter
%                                  (weight decay)
%                  'InputType'     [ {Direct} | SVDPreprocessed ] 
%                  'SingularValues'   { 3 } Number of singular values
%                                  maintained (in the case of SVD
%                                  preprocessing), ie. the number of
%                                  input neurons
%                  'Info'          [ {0} | 1 ] Running information
%                                  about the optimization process 
%
%       Output: E   Error
%               Y   Computed output
%               U   Network weights
%
%       Neural network regression is a generalization (with respect to
%       non-linearity) compared to the FIR model.
%
%       See also: lyngby_nn_qmain, lyngby_fir_main

% cvs : $Id: lyngby_nnr_main.m,v 1.8 1999/02/04 09:36:05 fnielsen Exp $
%       $Revision: 1.8 $


    % Default parameters
    validation = 1;
    sValidation = 'SingleBlocked';
    genOptim = 1;
    sGenOptim = 'Free';
    hiddenUnits = 3;
    Reg = 0.001;
    lag = 7;
    inputType = 1;
    singularValues = 3;
    bInfo = 1;
    spatialType = 1;
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));

      if strcmp(arg, 'genoptim')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  sGenOptim = arg;
	  if strcmp(arg, 'free')
	    genOptim = 1;
	  elseif strcmp(arg, 'earlystop')
	    genOptim = 2;
	  elseif strcmp(arg, 'hiddenunitsearlystop')
	    genOptim = 3;
	  elseif strcmp(arg, 'pruning')
	    genOptim = 4;
	  elseif strcmp(arg, 'pruning1dreggridsearch')
	    genOptim = 5;
	  elseif strcmp(arg, 'pruning2dreggridsearch')
	    genOptim = 6;
	  else
	    error(sprintf('Wrong arg to ''GenOptim'': %s', arg))
	  end
	else
	  error('The argument with ''GenOptim'' should be a string.'); 
	end

      elseif strcmp(arg, 'validation')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'singleblocked')
	    validation = 1;
	    sValidation = 'SingleBlocked';
	  else
	    error(sprintf('Wrong arg to ''Validation'': %s', arg))
	  end
	else
	  error('The argument with ''Validation'' should be a string.'); 
	end

      elseif strcmp(arg, 'hiddenunits')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  hiddenUnits = arg;
	else
	  error('The argument with ''HiddenUnits'' should be a number.'); 
	end
	
      elseif strcmp(arg, 'reg')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  Reg  = arg;
	else
	  error(['The argument with ''Reg'' should be a matrix with '  ...
		'1, 2 elements or a number corresponding to the number '...
		'of weights.']); 
	end
	
      elseif strcmp(arg, 'lag')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  lag  = arg;
	else
	  error('The argument with ''lag'' should be a number.'); 
	end
      
      elseif strcmp(arg, 'inputtype')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'direct')
	    inputType = 1;
	  elseif strcmp(arg, 'svdpreprocessed')
	    inputType = 2;
	  else
	    error('Argument with ''InputType'' is wrong');
	  end
	else
	  error('The argument with ''InputType'' should be a string.'); 
	end
      
      elseif strcmp(arg, 'singularvalues')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  singularValues  = arg;
	else
	  error('The argument with ''SingularValues'' should be a number.'); 
	end
      
      elseif strcmp(arg, 'info')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  bInfo = ~~arg;
	else
	  error('Argument to ''Info'' should be a value');
	end
	
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    % Neural network structure
    if inputType == 1
      Ni = lag;
    else
      Ni = singularValues;
    end
    Nii = Ni + 1;
    Nh = hiddenUnits;
    Nhh = Nh + 1;
    No = 1;
    Nu = Nii*Nh + Nhh*No;
    N = [ [Nii Nhh No] ; [Nii-1 Nh No] ];
    
    % Regularization
    % Reg = lyngby_nn_reg2reg(Reg, N);

    % Dimensions
    nVoxels = size(T,2);
    nSamples = size(x,1) - lag + 1;
    allSamples = size(x, 1);
    
    % Input
    Input = zeros(nSamples, lag);     % The Input as a matrix
    for r = 1:nSamples,
    	Input(r,:) = x(lag+r-1:-1:r)';
    end
    
    E = zeros(1,nVoxels);
    Y = zeros(nSamples, nVoxels);
    U = zeros( (lag+1)*hiddenUnits+hiddenUnits+1, nVoxels);
    
    if inputType == 2
      lyngby_log('SVD input to neural network');
      [Usvd, Ssvd, Vsvd] = svd(Input,0);
      Input = Usvd(:,1:singularValues);
      Vsvd = (Ssvd(1:singularValues, 1:singularValues) * ...
	    Vsvd(:, 1:singularValues)')';
    end
    
    if genOptim == 1
      TrainingIndices = 1:nSamples;
      ValIndices = [];
    elseif any(genOptim == [2 3 4 5 6])
      lyngby_log('Splitting dataset');
      nHalf = round(nSamples/2);
      TrainingIndices = 1:nHalf;
      ValIndices = (nHalf+1):nSamples;
    else
      error('Internal error')
    end
    
    
    if spatialType == 1
      % One single neural network for all voxels
      
      Target = T(lag-1+(1:nSamples),:);
      
      N = [ size(Input, 2)+1 hiddenUnits+1 1 ; ...
	size(Input, 2) hiddenUnits 1 ];
	
      WW = ones(1, nVoxels);
      [V, W] = lyngby_nn_initvw(N);
      
      maxSuperIterations = 10;
      for n = 1:maxSuperIterations
	
	TargetNN = Target / WW;
	normTargetNN = norm(TargetNN);
	TargetNN = TargetNN / normTargetNN;
	WW = WW * normTargetNN;
      
	[Eall, Yall, V, W] = lyngby_nn_qtrain(Input, TargetNN, ...
	    V, W, Reg, ...
	    'MaxIteration', 30, ...
	    'Info', bInfo);
	Yall = lyngby_nn_qforward(Input, V, W);
	WW = Yall \ Target;

	lyngby_log(sprintf('Neural training: %d / %d', n, ...
	    maxSuperIterations))

      end

      for v = 1:nVoxels
	WWW = W * WW(v);
	U(:,v) = [V(:) ; WWW(:)];
	Y(:,v) = lyngby_nn_qforward(Input, V, WWW);
	E(1,v) = lyngby_nn_qerror(Target(:,v), Y(:,v));
      end

    else
      % individual neural networks for each voxel
      
      for v = 1:nVoxels
	
	target = T(lag-1+(1:nSamples),v);
	
	[V, W] = lyngby_nn_qmain(Input, target, ...
	    'GenOptim', sGenOptim, ...
	    'Validation', sValidation, ...
	    'HiddenUnits', hiddenUnits, ...
	    'Reg', Reg, ...
	    'TrainSet', TrainingIndices, ...
	    'ValSet', ValIndices, ...
	    'Info', bInfo);
	Output = lyngby_nn_qforward(Input, V, W);
	E(1,v) = lyngby_nn_qerror(target, Output);
	
	if inputType == 2
	  V(1:Ni,:) = Vsvd * V(1:Ni,:)
	end
	
	U(:,v) = [V(:) ; W(:)];
	Y(:,v) = Output;
	
	%if ~lyngby_mod(v, 5)
	lyngby_log(sprintf('Neural training: %d / %d', v, ...
	    nVoxels))
	%end
	
      end
    end
      














