%FEATSELO Branch and bound feature selection
% 
% 	W = featselo(A,CRIT,K,T,FID)
%
% INPUT	
%   A     input dataset
%   CRIT  string name of the criterion or untrained mapping 
%           (optional, def= 'NN' 1-Nearest Neighbor error)
%   K     numner of features to select (optional, def: K=2)
%   T     validation set (optional)
%   N     Number of cross-validations (optional)
%   FID   File ID to write progress to (default [], see PRPROGRESS)
%
% OUTPUT
%   W     output feature selection mapping
% 
% DESCRIPTION
% Backward selection of K features by baktracking using the branch 
% and bound procedure on the data set A. CRIT sets the criterion 
% used by the feature evaluation routine FEATEVAL. If the data set T 
% is given, it is used as test set for FEATEVAL. Alternatively a number
% of cross-validations N may be supplied. The resulting W can be used for
% the selecting features of a dataset B by B*W. 
% 
% This procedure finds the optimum feature set if a monotoneous 
% criterion is used. The use of a testset does not guarantee that.
% 
% SEE ALSO
% MAPPINGS, DATASETS, FEATEVAL, FEATSELF, FEATSELB, FEATSELI,
% FEATSEL,  FEATSELP, FEATSELM, PRPROGRESS

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: featselo.m,v 1.10 2004/05/16 14:48:08 duin Exp $

function W = featselo(A,crit,kmin,T,fid)
	prtrace(mfilename);

	if nargin < 2 | isempty(crit), crit = 'NN'; end
	if nargin < 3 | isempty(kmin), kmin = 2; end
	if nargin < 4, T = []; end
	if (nargin < 5)
		fid = [];
	end

	if nargin == 0 | isempty(A)
      % Create an empty mapping:
		W = mapping('featselo',{crit,kmin,T});
		W = setname(W,'B&B FeatSel');
		return
	end

	isvaldset(A,1,2); % at least 1 object per class, 2 classes
	if ~isscalar(T), iscomdset(A,T); end

	[m,k,c] = getsize(A);
	featlist = getfeatlab(A);

	if ((kmin < 1) | (kmin >= k))
		error('The desired feature size should be > 0 and < dataset feature size')
	end
	
	% space for criteria values
	feat = zeros(1,k);

   % Get performance of the individual features:
	if isempty(T)
		for j=1:k
			feat(j) = feateval(A(:,j),crit);
		end
	elseif isscalar(T)
		for j=1:k
			feat(j) = feateval(A(:,j),crit,T);
		end
	else
		for j=1:k
			feat(j) = feateval(A(:,j),crit,T(:,j));
		end
	end

   % Get the kmin worst(?) individual features according to their
   % individual performance:
	[F,S] = sort(feat);
	Iopt = [k-kmin+1:k];

	I = [1:k];
	J = [zeros(1,kmin),1:(k-kmin-1),k-kmin+1,k+1];
	level = k;

   % Get the performance of Iopt
	if isdataset(T)
		bound = feateval(A(:,S(Iopt)),crit,T(:,S(Iopt)));
	elseif isscalar(T)
		bound = feateval(A(:,S(Iopt)),crit,T);
	else
		bound = feateval(A(:,S(Iopt)),crit);
	end

	C = inf;
	prprogress(fid,'\nfeatselo : Branch & Bound Feature Selection\n')
	while length(I) > 0 & J(k+1) == k+1;
		if J(level) == J(level+1) | level <= kmin | C <= bound
			J(level) = level - kmin;
			level = level + 1;
			I = sort([I,J(level)]);
			J(level) = J(level) + 1;
			C = inf;
		else
			I(J(level)) = [];
			level = level - 1;
			if J(level+1) < 3 & level == kmin+1 & 0 % never happens ??
				;
			else
				if isdataset(T)
					C = feateval(A(:,S(I)),crit,T(:,S(I)));
				elseif isscalar(T)
					C = feateval(A(:,S(I)),crit,T);
				else
					C = feateval(A(:,S(I)),crit);
				end
				prprogress(fid,'  level %i, crit: %10.5e',level,C);
				if level == kmin & C > bound
					bound = C;
					Iopt = I;
					prprogress(fid,'  FeatSet: ');
					prprogress(fid,' %i',S(I));
				end
				prprogress(fid,'\n')
			end
		end
	end
	prprogress(fid,'featselo  finished\n')

   % Store the optimal features in the mapping:
	W = featsel(k,S(Iopt));
	W = setlabels(W,featlist(S(Iopt),:));
	W = setname(W,'B&B FeatSel');

return
