source: prextra/lessc.m

function [w,I] = lessc(x, C, ftype, include_bias)
%LESSC Least Error in Sparse Subspaces classifier
%
%     W = LESSC(X, C, FTYPE, INCLUDE_BIAS)
%
% Train a linear classifier which also performs feature selection.
% In this version we do:
%                min \sum_i w_i + C*delta_i
%          s.t. forall_i   w^T f(x_i) > 1 - delta_i
%                       sum_i |w_i| = 1
% where f(x_i) is in principle free, but as a start we use the nearest
% mean idea:
%
%         f(x_i) = (x-mu2).^2 - (x-mu1).^2
% See for further definitions of f(x_i) lessfx.
%
% Dxd   15-3-2004
prtrace(mfilename);

if (nargin < 4)
        % To include a bias term in the model, we extend the number of features
        % by one:
        include_bias = 0;
end
if (nargin < 3)
        prwarning(3,'Use default function fx.');
        ftype = 1; 
end
if length(ftype)>1
        include_bias = ftype(2);
        ftype = ftype(1);
end
if (nargin < 2)
        prwarning(3,'C set to one');
        C = 1; 
end
if (nargin < 1) | (isempty(x))
        w = mapping(mfilename,{C,ftype,include_bias});
        w = setname(w,'LESS classifier');
        return
end

        
if ~ismapping(C)   % train the mapping

        % Unpack the dataset.
        islabtype(x,'crisp');
        isvaldset(x,1,2); % at least 1 object per class, 2 classes
        [m,k,c] = getsize(x); 

        if c == 2   % two-class classifier

                % get -1/+1 labels:
                nlab = getnlab(x);
                y = 2*nlab-3;

                % train and apply the feature mapping:
                par = lessfx(ftype,x);
                f = lessfx(par,x);

                if (include_bias)
                        f = [f ones(m,1)];
                        k = k+1;
                end

                % In the LP formulation, we define the free parameter vector as:
                %  [delta; w]
                % setup the constraints:
                yf = -repmat(y,1,k).*f;
                % standard version when we have Ax<b  and  Aeq x = b;
                A = [-eye(m)    -(+yf)];
                b = -ones(m,1);
                %Aeq = [zeros(1,m) ones(1,k)]; beq = 1;
                %if (include_bias), Aeq(1,end)=0; end
                Aeq = []; beq = [];
                % function to optimize:
                c = [repmat(C,1,m) ones(1,k)];
                %c = [ones(1,m) repmat(C,1,k)];
                if (include_bias), c(end) = 0; end
                % upper and lower bounds:
                lb = zeros(m+k,1);
                if (include_bias), lb(end) = -inf; end
                ub = repmat(inf,m+k,1);

                % optimize
                if (exist('glpkmex')==3)
                        [out,dummy]=glpkmex(1,c',A,b,repmat('U',m,1),lb,[],repmat('C',m+k,1));
                else
                        out = linprog(c,A,b,Aeq,beq,lb,ub);
                end
                w = out((m+1):end);

                % find out how many features are relevant:
                if (include_bias)
                        I = find(abs(w(1:(end-1)))>1e-8);
                        nr = length(I);
                else
                        I = find(abs(w)>1e-8);
                        nr = length(I);
                end

                % Store the classifier
                W.extend = include_bias;
                W.par = par;
                W.w = w;
                W.nr = nr;
                w = mapping(mfilename,'trained',W,getlablist(x),size(x,2),2);
                w = setname(w,'LESS classifier');
        
        else   % multi-class classifier:
                
                %error('Multiclass not implemented yet');
                w = mclassc(x,mapping(mfilename,{C,ftype,include_bias}));
                v = w.data{1}.data{1}.data.w;
                for i=2:length(w.data)
                        v = v + w.data{i}.data{1}.data.w;
                end
                I = find(abs(v)>0);
                
        end     
else
        % Evaluate the classifier on new data:
        W = getdata(C);

        % It is a simple linear classifier:
        if (W.extend)
                out = [lessfx(W.par,x) ones(size(x,1),1)]*W.w;
        else
                out = lessfx(W.par,x)*W.w;
        end

        % and put it nicely in a prtools dataset:
        w = setdat(x,sigm([out -out]),C);

end
                
return