function [V, W, EAcc, Info] = lyngby_nn_emain(X, T, arg1, arg2, arg3, ...
    arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, ...
    arg14, arg15, arg16)

% lyngby_nn_emain      - Main function for Entropic Neural network 
%
%	function [V, W] = lyngby_nn_emain(X, T, PropertyName, ...
%	    PropertyValue)
%
%       Input:    X   Neural network input
%                 T   Target output
%
%       Property: GenOptim      { {Free} | EarlyStop | 
%                               HiddenUnitsEarlyStop | Pruning 
%                               Pruning1DRegGridSearch |  
%                               Pruning2DRegGridSearch ]
%                               Generalization optimization
%                 Validation    [ {SingleBlocked | 
%                               BlockedCrossvalidation ]
%                               How to compute the generalization
%                               measure
%                 TrainSet      [ {not defined} | <indices> ]
%                               Indices for the training
%                               set. Default is all of X.
%                 ValSet        [ {not defined} | <indices> ]
%                               Indices for the validation
%                               set. Only used in connection with
%                               "non-free" generalization
%                               optimization. Default is no
%                               validation set.
%                 HiddenUnits   { 3 } Number of hidden units, not
%                               counting the threshold unit
%                 Reg           { 0.001 } Regularization parameter
%                               (weight decay)
%                 Seed          Seed for the random generator for
%                               the weights Matlab 4.x single seed
%                               generator is used. Default is no
%                               seed.
%                 Info          [ {0} | 1 ] Continuous information
%                               about the optimization
%
%       Output: V      Input weights
%               W      Output weights
%               EAcc   Evolution of error. Is dependent on the setting
%                      of 'GenOptim'
%               Info   Information about EAcc
%
%       Main function for the "entropic neural network" that is the
%       neural network with one binary output.
%
%       See also LYNGBY, LYNGBY_NN_QMAIN, LYNGBY_NN_CMAIN,
%                LYNGBY_NN_ETRAIN. 
%
% $Id: lyngby_nn_emain.m,v 1.24 2002/03/18 18:18:52 fnielsen Exp $


    % Check input
    if size(X,1) ~= size(T,1)
      error(sprintf(['''X'' and ''T'' should have the same number ' ...
	  'of examples (scans). X: %d, T: %d'], size(X,1), size(T,1)));
    elseif size(T, 2) ~= 1
      error(sprintf(['This neural network can only have on out: ' ...
	  '''T'' is too big. T: %d'], size(T,2))); 
    end

    % Default parameters
    genOptim = 1;          % Generalization optimization
    validation = 1;        % Generalization measure
    trainSet = 1:size(X,1);
    valSet = [];
    optMethod = 65;        % Optimization method
    
    Nh = 3;                % HiddenLayer
    [Np, Ni] = size(X);	   % Examples Total TrainingSet
    [Np, No] = size(T);
    Nii = Ni + 1;				

    Reg = 0.0001;          % Regularization, weight decay
    nPreEpochs	= 1300;	   % Pretrainingsepoch
    nReEpochs	= 300;	   % Retrainingepochs after pruning
    seed = []; 
    deleteFraction = 0.05;   % Pruning DeletePercent 
    nPruneEpochs = 120;    % Number of Prunings epochs
    bPruning     = 1;
    bInfo        = 0;
    minGradient  = 1e-6;
    maxIteration = 3000;
    
    % 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)
	  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, 'optmethod')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isnumeric(arg)
	  optMethod = arg;
	else
	  error('The argument with ''OptMethod'' should be a interger.'); 
	end

      elseif strcmp(arg, 'validation')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'singleblocked')
	    validation = 1;
	  elseif strcmp(arg, 'blockedcrossvalidation')
	    validation = 2;
	  else
	    error(sprintf('Wrong arg to ''Validation'': %s', arg))
	  end
	else
	  error('The argument with ''Validation'' should be a string.'); 
	end
	
      elseif strcmp(arg, 'trainset')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if max(arg) <= size(X,1) & min(arg) > 0
	    trainSet  = arg;
	  else
	    error(sprintf(['The indices for ''TrainSet'' should ' ...
		'be larger than 0 and not larger than X. ' ...
		'Max index: %d, X: %d'], max(arg), size(X,1)));
	  end
	else
	  error(['The argument with ''TrainSet'' should be a ' ...
		'vector of indices.']); 
	end

      elseif strcmp(arg, 'valset')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if isempty(arg) | (max(arg) <= size(X,1) & min(arg) > 0)
	    valSet  = arg;
	  else
	    error(sprintf(['The indices for ''valSet'' should ' ...
		'be larger than 0 and not larger than X. ' ...
		'Max index: %d, X: %d'], max(arg), size(X,1)));
	  end
	else
	  error(['The argument with ''valSet'' should be a ' ...
		' vector containing indices.']); 
	end

      elseif strcmp(arg, 'saveweights')
	% Njah, what should we do about the weights
	% Return 'em all, - that will be big.
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isstr(arg)
	  savevwFilename = arg;
	else
	  error('Argument with ''SaveWeights'' should be a string (a filename)');
	end
	
      elseif strcmp(arg, 'hiddenunits')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  Nh  = 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, 'info')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  bInfo  = ~~arg;
	else
	  error('The argument to ''Info'' should be a number.'); 
	end

      
      elseif strcmp(arg, 'seed')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  seed  = arg;
	else
	  error('The argument with ''Seed'' should be a number.'); 
	end

      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    % Sizes
    Nhh = Nh + 1;
    Nv  = Nii * Nh;
    Nw  = Nhh * No;
    Nu  = Nv + Nw;
    
    X	= [X ones(Np, 1)];
    
    % Weights 
    N = [ [Nii Nhh No] ; [Nii-1 Nh No] ];
    [V, W] = lyngby_nn_initvw(N, 'Seed', seed, 'Input', X);
    VMask = ones(size(V));
    WMask = ones(size(W));

    % Regularization
    if ~any(genOptim == [3 5 6])
      Reg = lyngby_nn_reg2reg(Reg, N);
    end

    % Construct training and validation set
    if genOptim ~= 1
      if isempty(valSet)
	% Contruct validation set
	if isempty(trainSet)
	  % Divide
	  error(['Not implemented yet: Automatic train and ' ...
		'validation set']);
	else
	  % Use the rest for validation set.
	  valSet = lyngby_set_diff(1:size(X,1), trainSet); 
	end
      else
	if isempty(trainSet)
	  % Use the rest for training set.
	  trainSet = lyngby_set_diff(1:size(X,1), valSet); 
	elseif ~isempty(lyngby_set_and(valSet,trainSet))
	  error(sprintf(['Training set and Validation set should ' ...
		'not be overlaping.']));
	end
      end
    end

    
    % Training and Validation Set
    X0 = X(trainSet,:);
    T0 = T(trainSet,:);
    if ~isempty(valSet)
      X1 = X(valSet,:);
      T1 = T(valSet,:);
    end

    % Generalization/network Optimization
    if genOptim == 1
      % Free
      [E, Y, V, W] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	  'MaxIteration', maxIteration, ...
	  'MinGradient', minGradient, ...
	  'Method', optMethod, ...
	  'Info', bInfo);
      if ~isempty(valSet)
	% If E should be returned...
	[Y1, H] = lyngby_nn_eforward(X1, V, W);
	E      = lyngby_nn_eerror(T1, Y1);
      end

      EAcc = E;
      Info = [];
      
    elseif genOptim == 2
      % Early stop

      cont = 1;
      while cont
	[E0, Y0, VAcc, WAcc] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	    'MaxIteration', maxIteration, ...
	    'MinGradient', minGradient, ...
	    'WeightAcc', 1, ...
	    'Method', optMethod, ...
	    'Info', bInfo);
      
	lyngby_log('Computing validation error');
	
	N = [[Nii Nhh 0] ; [ Ni Nh No]];
	E1Acc = zeros(size(VAcc,2), 1);
	for e = 1:size(VAcc,2)
	  [V, W]   = lyngby_nn_u2vw([VAcc(:,e) ; WAcc(:,e)], N);
	  [Y1, H]  = lyngby_nn_eforward(X1, V, W);
	  EAcc(e) = lyngby_nn_eerror(T1, Y1);
	end

	% Find lowest validation error
	[E, Eindex] = min(EAcc);
	[V, W]   = lyngby_nn_u2vw([VAcc(:,Eindex) ; WAcc(:,Eindex)], N);
	
	if 1==1   % Something intelligent here !!!
	  cont = 0;
	end
      end
	
      Info = [];
	
    elseif genOptim == 3
      % HiddenUnitsEarlyStop

      RegOrig = Reg;
      maxNh = Nh;
      
      NhArray = lyngby_set_unique([1 round(exp(mean(log([1 maxNh])))) ...
	    Nh]);

      EAcc = Inf * ones(MaxIteration, length(NhArray));
      
      Nh = 1;
      N = [ [Nii Nh+1 No] ; [Ni  Nh  No] ];
      [V, W] = lyngby_nn_initvw(N);
      Reg = lyngby_nn_reg2reg(RegOrig, N);
      lyngby_log(sprintf('Optimizing for %d hidden units', Nh)); 
      [E0, Y0, VAcc1, WAcc1] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	  'MaxIteration', maxIteration, ...
	  'MinGradient', minGradient, ...
	  'WeightAcc', 1, ...
	  'Method', optMethod, ...
	  'Info', bInfo);
      
      lyngby_log('Computing validation error');
      for e = 1:size(VAcc1,2)
	[V, W]   = lyngby_nn_u2vw([VAcc1(:,e) ; WAcc1(:,e)], N);
	[Y1, H]  = lyngby_nn_eforward(X1, V, W);
	EAcc(e, 1) = lyngby_nn_eerror(T1, Y1);
      end

      % Find lowest validation error
      [E, indexEpoch] = min(EAcc);
      
      if length(NhArray) == 1
	[V, W] = lyngby_nn_u2vw([VAcc(:, indexEpoch) ; WAcc(:,indexEpoch)], N);
      else
	Nh = maxNh;
	N = [ [Nii Nh+1 No] ; [Nii-1 Nh No] ];
	[V, W] = lyngby_nn_initvw(N);
	Reg = lyngby_nn_reg2reg(RegOrig, N);
	lyngby_log(sprintf('Optimizing for %d hidden units', Nh)); 
	[E0, Y0, VAcc2, WAcc2] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	    'MaxIteration', maxIteration, ...
	    'MinGradient', minGradient', ...
	    'WeightAcc', 1, ...
	    'Method', optMethod, ...
	    'Info', bInfo);
      
	lyngby_log('Computing validation error');
	for e = 1:size(VAcc2,2)
	  [V, W]   = lyngby_nn_u2vw([VAcc2(:,e) ; WAcc2(:,e)], N);
	  [Y1, H]  = lyngby_nn_eforward(X1, V, W);
	  EAcc(e, 2) = lyngby_nn_eerror(T1, Y1);
	end
	
      end
      if length(NhArray) == 2 
	[E, indexEpoch] = min(EAcc);
	[E, indexNh] = min(E);

	if indexNh == 1
	  Nh = 1;
	  N = [ [Nii Nh+1 No] ; [Nii-1 Nh No] ];
	  [V, W]   = lyngby_nn_u2vw([VAcc1(:,indexEpoch(1)) ; ...
	    WAcc1(:,indexEpoch(1))], N);
	else
	  [V, W]   = lyngby_nn_u2vw([VAcc2(:,indexEpoch(2)) ; ...
	    WAcc2(:,indexEpoch(2))], N);
	end
	  
      else

	Nh = NhArray(2);
	cont = 1;
	while cont
	  N = [ [Nii Nh+1 No] ; [Nii-1  Nh  No] ];
	  [V, W] = lyngby_nn_initvw(N);
	  Reg = lyngby_nn_reg2reg(RegOrig, N);
	  lyngby_log(sprintf('Optimizing for %d hidden units', Nh)); 
	  [E0, Y0, VAcc3, WAcc3] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	      'MaxIteration', maxIteration, ...
	      'MinGradient', minGradient, ...
	      'WeightAcc', 1, ...
	      'Method', optMethod, ...
	      'Info', bInfo);
	  
	  lyngby_log('Computing validation error');
	  for e = 1:size(VAcc3,2)
	    [V, W]   = lyngby_nn_u2vw([VAcc3(:,e) ; WAcc3(:,e)], N);
	    [Y1, H]  = lyngby_nn_eforward(X1, V, W);
	    EAcc(e, 3) = lyngby_nn_eerror(T1, Y1);
	  end
	  [E, indexEpoch] = min(EAcc);
	  [E, indexNh] = min(E);

	  if indexNh == 1
	    Nh = 1;
	    N = [ [Nii Nh+1 No] ; [Nii-1 Nh No] ];
	    [V, W]   = lyngby_nn_u2vw([VAcc1(:,indexEpoch(1)) ; ...
	      WAcc1(:,indexEpoch(1))], N);
	  elseif indexNh == 2
	    Nh = maxNh;
	    N = [ [Nii Nh+1 No] ; [Nii-1 Nh No] ];
	    [V, W]   = lyngby_nn_u2vw([VAcc2(:,indexEpoch(2)) ; ...
	      WAcc2(:,indexEpoch(2))], N);
	  else
	    [V, W]   = lyngby_nn_u2vw([VAcc3(:,indexEpoch(3)) ; ...
	      WAcc3(:,indexEpoch(3))], N);
	  end
	  
	  cont = 0;
	  
	end
      end
	
      [Info, index] = sort(-NhArray);
      Info = -Info;
      EAcc = EAcc(:,index);
      
    elseif genOptim == 4
      % Pruning
      
      % Pretraining
      lyngby_log('Pretraining before pruning')
      [E,Y,V,W] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	    'MaxIteration', maxIteration, ...
	    'MinGradient', minGradient, ...
	    'Method', optMethod, ...
	    'Info', bInfo);
      
      % Pruning loop
      e = 1;
      while 1==1
	
	lyngby_log('Computing validation error');
	[Y1, H]  = lyngby_nn_eforward(X1, V, W);
	EAcc(e,1) = lyngby_nn_eerror(T1, Y1);
	Info(e,1) = sum(V(:) ~= 0) + sum(W(:) ~= 0);

	% Save weights
	[UAcc(:,e), N] = lyngby_nn_vw2u(V,W);
	
	lyngby_log('Pruning');
	[VMask,WMask] = lyngby_nn_epruneobd(X, T, V, W, Reg, ...
	    'DeleteFraction', deleteFraction);
	V = V .* VMask;
	W = W .* WMask;
	
	% Collaps V, W, X and Reg
	
	if (sum(sum(VMask)) + sum(sum(WMask))<=1)
	  break
	end

	lyngby_log(sprintf('Retraining for %d parameters', ...
	    sum([VMask(:) ; WMask(:)])));
	[E,Y,V,W] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	    'MaxIteration', maxIteration, ...
	    'MinGradient', minGradient, ...
	    'Method', optMethod, ...
	    'Info', bInfo);
	
	e = e+1;
      end
      
      % Return best V Wfrom saved weights
      [dummy, index] = min(EAcc);
      [V,W] = lyngby_nn_u2vw(UAcc(:,index), N);
      
    elseif genOptim == 5
      % 1D grid search for regularization parameter and pruning

      if length(Reg) ~= 2
	error(['In the case of ''Pruning / 1D reg. grid search'' ' ...
	      'the regularization should be of length 2']);
      elseif Reg(1) >= Reg(2)
	error(['The first element in the regularization should be ' ...
	      'larger than the second.']);
      end
      
      % Regularization
      steps = 10;
      lReg = log(Reg);
      regV = exp(lReg(1) : (lReg(2)-lReg(1))/(steps-1) : lReg(2));
      regW = regV;
      
      Info = zeros(steps,3);
      EAcc = zeros(steps,1);
      
      V0 = V;
      W0 = W;
      for r = 1:steps
	  
	% Weights
	V = V0;
	W = W0;
	
	% Regularization 
	Reg = [regV(r) regW(r)];
	Info(r,1:2) = Reg; 
	Reg = lyngby_nn_reg2reg(Reg, N);
	
	% Pretraining
	lyngby_log(sprintf('Pretraining: %d/%d', r, steps));
	[E,Y,V,W] = lyngby_nn_etrain(X0, T0, V0, W0, Reg, ...
	    'MaxIteration', maxIteration, ...
	    'MinGradient', minGradient, ...
	    'Method', optMethod, ...
	    'Info', bInfo);      
	
	% Pruning loop
	EAccPrune = [];
	UAccPrune = [];
	e = 1;
	while 1==1
	  
	  lyngby_log('Computing validation error');
	  [Y1, H]  = lyngby_nn_eforward(X1, V, W);
	  EAccPrune(e,1) = lyngby_nn_eerror(T1, Y1);
	  
	  % Save weights
	  [UAccPrune(:,e), N] = lyngby_nn_vw2u(V,W);
	  
	  lyngby_log('Pruning');
	  [VMask,WMask] = lyngby_nn_epruneobd(X, T, V, W, Reg, ...
	      'DeleteFraction', deleteFraction);
	  V = V .* VMask;
	  W = W .* WMask;
	  
	  % Collaps V, W, X and Reg
	  
	  if (sum(sum(VMask)) + sum(sum(WMask))<=1)
	    break
	  end
	  
	  lyngby_log(sprintf('Retraining: %d/%d for %d parameters', ...
	      r, steps, sum([VMask(:) ; WMask(:)])));
	  [E,Y,V,W] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
	      'MaxIteration', maxIteration, ...
	      'MinGradient', minGradient, ...
	      'Method', optMethod, ...
	      'Info', bInfo);
	  
	  e = e+1;
	end
	
	% Return best pruned network with a specific set of
	% regarization parameters
	[EAcc(r), index] = min(EAccPrune);
	UAcc(:,r) = UAccPrune(:,index);
	
      end

      % Return best overall network
      [dummy, index] = min(EAcc);
      [V,W] = lyngby_nn_u2vw(UAcc(:,index), N);
    
    elseif genOptim == 6
      % 2D grid search for regularization parameter and pruning

      if length(Reg) ~= 4
	error(['In the case of ''Pruning / 2D reg. grid search'' ' ...
	      'the regularization should be of length 4']);
      elseif Reg(1) >= Reg(2)
	error(['The first element in the regularization should be ' ...
	      'larger than the second.']);
      elseif Reg(3) >= Reg(4)
	error(['The third element in the regularization should be ' ...
	      'larger than the second.']);
      end
      
      % Regularization
      steps = 10;
      lReg = log(Reg);
      regV = exp(lReg(1) : (lReg(2)-lReg(1))/(steps-1) : lReg(2));
      regW = exp(lReg(3) : (lReg(4)-lReg(3))/(steps-1) : lReg(4));
      

      Info = zeros(steps*steps,3);
      EAcc = zeros(steps*steps,1);
      
      V0 = V;
      W0 = W;
      for r1 = 1:steps
	for r2 = 1:steps
	  r = steps*(r1-1)+r2;
	  
	  % Weights
	  V = V0;
	  W = W0;
	  
	  % Regularization 
	  Reg = [regV(r1) regW(r2)];
	  Info(r,1:2) = Reg; 
	  Reg = lyngby_nn_reg2reg(Reg, N);
	  
	  % Pretraining
	  lyngby_log(sprintf('Pretraining: %d/%d-%d/%d', ...
	      r1, steps, r2, steps));
	  [E,Y,V,W] = lyngby_nn_etrain(X0, T0, V0, W0, Reg, ...
	      'MaxIteration', maxIteration, ...
	      'MinGradient', minGradient, ...
	      'Method', optMethod, ...
	      'Info', bInfo);      
	  
	  % Pruning loop
	  EAccPrune = [];
	  UAccPrune = [];
	  e = 1;
	  while 1==1
	    
	    lyngby_log('Computing validation error');
	    [Y1, H]  = lyngby_nn_eforward(X1, V, W);
	    EAccPrune(e,1) = lyngby_nn_eerror(T1, Y1);
	    
	    % Save weights
	    [UAccPrune(:,e), N] = lyngby_nn_vw2u(V,W);
	    
	    lyngby_log('Pruning');
	    [VMask,WMask] = lyngby_nn_epruneobd(X, T, V, W, Reg, ...
		'DeleteFraction', deleteFraction);
	    V = V .* VMask;
	    W = W .* WMask;
	    
	    % Collaps V, W, X and Reg
	    
	    if (sum(sum(VMask)) + sum(sum(WMask))<=1)
	      break
	    end

	    lyngby_log(sprintf('Retraining: %d/%d-%d/%d for %d parameters', ...
		r1, steps, r2, steps, sum([VMask(:) ; WMask(:)])));
	    [E,Y,V,W] = lyngby_nn_etrain(X0, T0, V, W, Reg, ...
		'MaxIteration', maxIteration, ...
		'MinGradient', minGradient, ...
		'Method', optMethod, ...
		'Info', bInfo);
	
	    e = e+1;
	  end
      
	  % Return best pruned network with a specific set of
	  % regarization parameters
	  [EAcc(r), index] = min(EAccPrune);
	  UAcc(:,r) = UAccPrune(:,index);
	  
	end
      end

      % Return best overall network
      [dummy, index] = min(EAcc);
      [V,W] = lyngby_nn_u2vw(UAcc(:,index), N);
    
    end










