function [VMask,WMask] = lyngby_nn_pruneobd(V, W, ddCV, ddCW, arg1, ...
    arg2)

% lyngby_nn_pruneobd   - Pruning by Optimal Brain Damage
%
%	function [VMask,WMask] = lyngby_nn_pruneobd(V, W, ddCV,ddCW,
%           'PropertyName', 'PropertyValue')
%
%	Input:	V      Input weights
%               W      Output weights
%               ddCV   Second order derivative of the costfunction with
%                      respect to the input weights, matrix
%                      conataining the diagonal of the hessian
%               ddCW   Second order derivative of the costfunction
%                      with respect to the output weights, matrix
%                      containing the diagonal of the hessian
%               Property:
%                  'DeleteFraction'   {0.05} Fraction to delete 
%                                     of remaining weights
%
%	Output:	VMask   Mask for input weights. Deleted weight is
%                       zero, else one
%               WMask   Mask for output weights.
%
%       Optimal Brain Damage pruning for a neural network.
%
%       See also: lyngby_nn_cpruneobd, lyngby_nn_eprunobd,
%                 lyngby_nn_qpruneobd

% cvs : $Id: lyngby_nn_pruneobd.m,v 1.4 1998/04/16 11:15:57 fnielsen Exp $
%       $Revision: 1.4 $

    % Default
    deleteFraction = 0.05;

    % Parse Properties
    n = 1;
    while n <= nargin-6
      eval(sprintf('arg = lower(arg%d);', n));

      if strcmp(arg, 'DeleteFraction')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isreal(arg)
	  if arg > 0 & arg < 1
	    deleteFraction = arg;
	  else
	    error(sprintf('Wrong arg to ''DeleteFraction'': %s', arg))
	  end
	else
	  error('The argument with ''DeleteFraction'' should be a string.'); 
	end
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    VMask = V~=0;
    WMask = W~=0;
    
    [Nii, Nh] = size(V);
    [Nhh, No] = size(W);
    Ni = Nii - 1;
    Nv = Nii * Nh;
    Nw = Nhh * No;
    Nu = Nv + Nw;

    % Saliency
    SalV = ddCV .* V.^2; 
    SalW = ddCW .* W.^2;
    
    % Find the number of weights to delete
    NBack   = prod(size(find(V))) + prod(size(find(W)));
    NDeleted = Nu - NBack;
    NDelete  = ceil(NBack*deleteFraction);
    
    % Nonzero weights with zero saliency is a problem in connection
    % with sorting
    SalV(find(V~=0 & SalV==0)) = realmin;
    SalW(find(W~=0 & SalW==0)) = realmin;
    
    % Find weights with minimum saliency	
    SalList = [ SalV(:) ; SalW(:) ];
    [SalListSort,SalListIndex] 	= sort(SalList);	

    for u = (1+NDeleted) : (NDeleted+NDelete)
      if SalListIndex(u) <= Nv
	% Delete in InputHiddenLayer (zero a weight)
	h = fix((SalListIndex(u)-1)/Nii) + 1;
	i = SalListIndex(u)-(h-1)*Nii;
	VMask(i, h) = 0;
	if ~any(VMask(:,h))
	  WMask(h,:) = 0;
	end
      else
	% Delete in HiddenOutputLayer 
	o = fix((SalListIndex(u)-Nv-1) / Nhh) + 1;	% The to-be deleted unit
	h = (SalListIndex(u)-Nv) - (o-1)*Nhh;
	WMask(h, o) = 0;
	if h ~= Nhh & ~any(WMask(h,:))
	  VMask(:,h) = 0;
	end
      end
    end











