function [S, S1, Y, V, W, EAcc, Info] = lyngby_nns_main(X, T, arg1, arg2, ...
    arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, ...
    arg13, arg14, arg15, arg16, arg17, arg18, arg19, arg20, arg21, arg22)

% lyngby_nns_main      - Main function for Neural network saliency
%
%	function [S, S1, Y, V, W, EAcc, Info] = lyngby_nns_main(X, T,
%	    'PropertyName', 'PropertyValue')
%
%       Input:  X   The datamatrix
%               T   Target output (the paradigm)
%               PropertyName:
%                  'Type'          [ {LS} | LSSecondOrder | 
%                                  Entropy | EntropySecondOrder ]
%                  'SVDComponents'   { 50 } Number of singular values
%                                  maintained, ie. the number of
%                                  neural network inputs.
%                  'NNInput'       [ {SVD} | SVDNormalized | Direct |
%                                  SOP | SOPNormalized | SOPResSVD ]
%                                  Type of signal on the neural
%                                  network input
%                  'HiddenUnits'   { 3 } Number of hidden units, not
%                                  counting the threshold unit
%                  'Reg'           { 0.001 } Regularization parameter
%                                   (weight decay)
%                  'GenOptim'      [ {Free} | EarlyStop |
%                                  HiddenUnitsEarlyStop | Pruning |
%                                  Pruning1DRegGridSearch |
%                                  Pruning2DReggridSearch ]
%                                  Generalization optimization
%                  'Validation'    [ {InBasisSingleBlocked} | 
%                                  OffBasisSingleBlocked ]
%                                  Type of validation (the type of
%                                  validation set)
%                  'Run'           Run specification, used in the
%                                  automatic splitting of training and
%                                  validation set. 
%                  'Info'          [ {0} | 1 ] Continuous information
%                                  about the optimization
%
%       Output: S      Saliency
%               S1     First order saliency
%               E      Error for each example (each scan)
%               V      Input weights
%               W      Output weights
%               EAcc   Is the accumulated error (the evolution of the
%                      error), and it will depend on the setting of
%                      'GenOptim'
%               Info   Information about EAcc
%
%       This function will first SVD the 'X' matrix into 'SVDComponents'
%       singular values and rotations. The output from the SVD will be
%       feed to a two-layer feed-forward neural network. This neural
%       network will be optimized for the most generalizing prediction
%       of 'T' (the paradigm) according to the setting of 'GenOptim'
%       and 'Validation'. After this optimization the saliency of the
%       variables in 'X' will be calculated.
%
%       See also: lyngby_svd, lyngby_nn_emain, lyngby_nn_qmain

% cvs : $Id: lyngby_nns_main.m,v 1.21 1998/09/22 14:05:54 fnielsen Exp $
%       $Revision: 1.21 $

    lyngby_global

    % Sizes
    [rX, cX] = size(X);
    [rT, cT] = size(T);
    
    if rX ~= rT
      error(sprintf(['X and T do not have the same length. X: %d, ' ...
	'T: %d'], rX, rT));
    end
    
    % Default parameters
    saliencyType = 1;
    svdComponents = 50;
    canComponents = 10;
    nnInput = 1;    
    validation = 1;
    sValidation = 'SingleBlocked';
    genOptim = 1;
    sGenOptim = 'Free';
    hiddenUnits = 3;
    Reg = 0.001;
    Run = [];
    bInfo = 0;
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'type')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'ls')
	    saliencyType = 1;
	  elseif strcmp(arg, 'lssecondorder')
	    saliencyType = 2;
	  elseif strcmp(arg, 'entropy')
	    saliencyType = 3;
	  elseif strcmp(arg, 'entropysecondorder')
	    saliencyType = 4;
	  else
	    error(sprintf('Wrong arg to ''Type'': %s', arg))
	  end
	else
	  error('The argument with ''Type'' should be a string.'); 
	end
      elseif strcmp(arg, 'svdcomponents')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  svdComponents  = arg;
	else
	  error('The argument with ''SVDComponents'' should be a number.'); 
	end
	
      elseif strcmp(arg, 'nninput')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'svd')
	    nnInput = 1;
	  elseif strcmp(arg, 'svdnormalized')
	    nnInput = 2;
	  elseif strcmp(arg, 'direct')
	    nnInput = 3;
	  elseif strcmp(arg, 'sop')
	    nnInput = 4;
	  elseif strcmp(arg, 'sopnormalized')
	    nnInput = 5;
	  elseif strcmp(arg, 'sopressvd')
	    nnInput = 6;
	  else
	    error(sprintf('Wrong arg to ''NNInput'': %s', arg))
	  end
	else
	  error('The argument with ''NNInput'' 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, 'inbasissingleblocked')
	    validation = 1;
	    sValidation = 'SingleBlocked';
	  elseif strcmp(arg, 'offbasissingleblocked')
	    validation = 2;
	    sValidation = 'SingleBlocked';
	  %elseif strcmp(arg, 'inbasisblockedcrossvalidation')
	  %  validation = 3;
	  %  sValidation = 'BlockedCrossValidation';
	  else
	    error(sprintf('Wrong arg to ''Validation'': %s', arg))
	  end
	else
	  error('The argument with ''Validation'' should be a string.'); 
	end
	
      elseif 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, '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, 'run')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if size(arg, 1) == rX
	    Run = arg;
	  else
	    error(sprintf(['The argument with ''Run'' should have ' ...
		  'the same number of rows as X, ', ...
		  'X: %dx%d, Run: %dx%d'], size(X), size(arg)));
	  end
	else
	  error(['The argument with ''Run'' should contain numbers']);
	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

    if any(nnInput == [ 4 5 6 ])
      if isempty(Run)
	error(['When the input to the neural network is SOP ', ...
	      'then the ''Run'' property should be defined']);
      end
    end
    
    Nl = size(X, 2);
    Ni = svdComponents;
    Nii = Ni + 1;
    Nh = hiddenUnits;
    Nhh = Nh + 1;
    No = size(T,2);
    Np = size(X,1);

    lyngby_log('Initializing ...');
    
    nVoxels = size(X,2);
    nSamples = size(X,1);

    % Split training and test set
    if genOptim == 1
      TrainIndices = 1:size(X,1);
      ValIndices = [];
    else
      [TrainIndices, ValIndices] = lyngby_nn_setsplit( ...
	  'Run', Run, 'Paradigm', T);
    end

    % Projection 
    if genOptim == 1 | any(validation == [1 2])
      % In basis validation
      
      lyngby_log('Computing SVD ...');
      
      if nnInput == 3
	% Direct
	I = X;
	Vsvd = eye(size(X, 2));
      elseif any(nnInput==[1 2])
	% SVD
	if validation == 1
	  [Usvd, Ssvd, Vsvd] = lyngby_svd(X, 'ReduceComponents', ...
	      svdComponents);
	else
	  [Usvd, Ssvd, Vsvd] = lyngby_svd(X(TrainIndices,:), ...
	      'ReduceComponents', svdComponents);
	end
	if nnInput == 1
	  % SVD direct
	  I = Usvd * diag(Ssvd);
	elseif nnInput == 2
	  % SVD normalized
	  I = Usvd;
	  Vsvd =  Vsvd * diag(1./Ssvd);
	else
	  error('Ops, wrong nninput');
	end
	% Project validation set
	if validation == 2
	  I = [ I ; X(ValIndices,:) * Vsvd ];
	end
      elseif any(nnInput==[4 5])
	% SOP
	if validation == 1
	  [Usvd, Ssvd, Vsvd] = lyngby_sop_main(X, Run, ...
	      'Components', svdComponents);
	else
	  [Usvd, Ssvd, Vsvd] = lyngby_sop_main(X(TrainIndices,:), ...
	      Run(TrainIndices,:), 'Components', svdComponents);
	end
	if nnInput == 4
	  % SOP direct 
	  I = X * Vsvd;
	elseif nnInput == 5  
	  % SOP normalized
	  Vsvd =  Vsvd * diag(1./Ssvd);
	  I = X * Vsvd;
	end
      elseif nnInput == 6
	% SOPResSVD
	
	[Usvd,Ssvd,Vsvd] = lyngby_sopressvd(X(TrainIndices,:), ...
	    Run(TrainIndices,:), ...
	    'CanComponents', canComponents, ...
	    'PrincComponents', svdComponents);
	I = X * Vsvd;
	
      else
	error('Ops, wrong nninput');
      end
    else
      error('Ops, no projection ');
    end
    
    if any(saliencyType == [1 2])
      % ls, Quadratic
      
      lyngby_log('Optimizing neural network (quadratic) ...');

      [V, W, EAcc, Info] = lyngby_nn_qmain(I, T, ...
	  'GenOptim', sGenOptim, ...
	  'Validation', sValidation, ...
	  'TrainSet', TrainIndices, 'ValSet', ValIndices, ...
	  'Reg', Reg, ...
	  'HiddenUnits', hiddenUnits, ...
	  'Info', bInfo);
      
      % Nh might be changed:
      [Nii, Nh] = size(V);
      Ni = Nii - 1;
      
      [Y, H] = lyngby_nn_qforward(I, V, W);
      E = lyngby_nn_qerror(T, Y);
      
      lyngby_log('Computing saliency (quadratic) ...');

      if saliencyType == 1
	% saliency, quadratic, first order, type 1
	
	S1 = 1/(No*Np) * ...
	    sum(X .* ...
	    (( ((T-Y) * W(1:Nh,:)') .* (1-H(:,1:Nh).^2) ) * ...
	    (V(1:Ni,:)' * Vsvd')));
	S = S1;
	
      else
	% saliency, quadratic, first and second order, type 1
	S1 = 1/(No*Np) * ...
	    sum(X .* ...
	    (( ((T-Y) * W(1:Nh,:)') .* (1-H(:,1:Nh).^2) ) * ...
	    (V(1:Ni,:)' * Vsvd')));
	if No == 1
	  S2 = 1/(2*No*Np) * ...
	      sum(X.^2 .* ...
	      (((1-H(:,1:Nh).^2) * ...
	      (W(1:Nh,ones(1,Ni)) .* V(1:Ni,:)')).^2 * ...
	      Vsvd'.^2));
	else
	  indexO = kron(ones(1,Ni*Np), 1:No);
	  indexP = kron(ones(1,Ni), kron(1:Np, ones(1,No)));
	  indexI = kron(1:Ni, ones(1,No*Np));
	  indexH = 1:Nh;
	  S2 = 1/(2*No*Np) * ...
	      sum(X.^2 .* ...
	      (reshape( ...
	      sum(reshape(sum( W(indexH,indexO) .* ...
	      (1-H(indexP,indexH)'.^2)  .* ...
	      V(indexI, indexH)').^2, No, Np*Ni)), ...
	      Np, Ni) * Vsvd'.^2));
	end

	S = S1 + S2;

	if 1==0
	  % The matrix form above should correspond to this:
	  for l = 1:Nl
	    S2test(l) = 0;
	    n = 1;
	    for p = 1:Np
	      for o = 1:No
		for i = 1:Ni
		  for h = 1:Nh
		    s1(h) = W(h,o) * (1-H(p,h)^2) * V(i,h); 
		  end
		  s2(n) = sum(s1)^2 * Vsvd(l,i)^2 * X(p,l)^2;
		  n = n+1;
		end
	      end
	    end
	    S2test(l) = 1/(2*No*Np) * sum(s2);
	    disp(sprintf('%d/%d', l, Nl));
	  end
	  % The difference should be in the area of 'eps'*S2
	  max(S2test-S2)
	end
      end
	
    elseif any(saliencyType == [3 4])
      T = lyngby_nn_etarget(T);
      
      lyngby_log('Optimizing neural network (entropic) ...');

      [V, W, EAcc, Info] = lyngby_nn_emain(I, T, ...
	  'GenOptim', sGenOptim, ...
	  'Validation', sValidation, ...
	  'TrainSet', TrainIndices, 'ValSet', ValIndices, ...
	  'Reg', Reg, ...
	  'HiddenUnits', hiddenUnits, ...
          'Info', bInfo);
      
      % Nh might be changed:
      [Nii, Nh] = size(V);
      Ni = Nii - 1;
      
      [Y, H] = lyngby_nn_eforward(I, V, W);
      E = lyngby_nn_eerror(T, Y);

      lyngby_log('Computing saliency (entropic) ...');

      if saliencyType == 3
	% saliency, entropy, first order, type 1
	
	S1 = 1/(No*Np) * ...
	    sum(X .* ...
	    (( ((T-Y) * W(1:Nh,:)') .* (1-H(:,1:Nh).^2) ) * ...
	    (V(1:Ni,:)' * Vsvd')));
	S = S1;
      
      elseif saliencyType == 4
	% saliency, entropy, first and second order, type 1

	S1 = 1/(No*Np) * ...
	    sum(X .* ...
	    (( ((T-Y) * W(1:Nh,:)') .* (1-H(:,1:Nh).^2) ) * ...
	    (V(1:Ni,:)' * Vsvd')));
	S2 = 1/(2*No*Np) * ...
	    sum(X.^2 .* ...
	    (((1-Y.^2) * ones(1,Ni)) .* ...
	    ((1-H(:,1:Nh).^2) * ...
	    (W(1:Nh,ones(1,Ni)) .* V(1:Ni,:)')).^2 * ...
	    Vsvd'.^2));
	S = S1 + S2;
      end
    else
      error('Internal Error');
    end
      





























