function model=ekozinec(data,options,init_model)% EKOZINEC Kozinec's algorithm for eps-optimal separating hyperplane.%% Synopsis:%  model = ekozinec(data)%  model = ekozinec(data,options)%  model = ekozinec(data,options,init_model)%% Description:%  This function is implementation of the Kozinec's algorithm%  with eps-optimality stopping condition [SH10]. The algorithm %  finds the eps-optimal separating hyperplane.% %  model=ekozinec(data) the Kozinec's rule is used to find the closest %   points w1, w2 from the convex hulls of the vectors from the first and %   the second class. The found points determine the optimal separating %   hyperplane. % %  model=ekozinec(data,options) specifies stopping conditions of%   the algorithm in structure options:%    .eps [1x1] ... controls how close is the found solution to%     the optimal hyperplane in terms of margin %     (default eps=0.01). The options for eps are: %       eps > 0 ... eps-optimal hyperplane is sought for.%       eps == 0 ... algorithm converges to the optimal hyperplane (but it%                    does not have to stop in finite number of iterations).%       eps < 0 ... algorithm stops when the separating hyperplane %                   is found (zero training error) regardless the margin %                   so it solves the same task as the ordinary Perceptron.%    .tmax [1x1]... maximal number of iterations.%%  model = ekozinec(data,options,init_model) specifies initial model%   which must contain:%    .W1 [dim x 1] ... Vector from the first convex hull.%    .W2 [dim x 1] ... Vector from the second convex hull.%% Input:%  data [struct] Labeled (binary) training data. %   .X [dim x num_data] Input vectors.%   .y [1 x num_data] Labels (1 or 2).%%  options [struct] %   .eps [real] Controls how closeness to the optimal hypeprlane (see above).%   .tmax [1x1] Maximal number of iterations (default tmax=inf).%  %  init_model [struct] Initial model; must contain items%    .W1 [dim x 1], .W2 [dim x 1] see above.%% Output:%  model [struct] Binary linear classifier:%   .W [dim x 1] Normal vector of hyperplane.%   .b [1x1] Bias of hyperplane.%  %   .W1 [dim x 1] The nearest vector of the first convex hull.%   .W2 [dim x 1] The nearest vector of the second convex hull.%   .margin [1x1] Margin of the found hyperplane.%   .exitflag [1x1] 1 ... eps-optimality condition satisfied or separating%                         hyperplane has been found %                   0 ... number of iterations exceeded tmax.%   .t [1x1] Number of iterations.%% Example:%  data = genlsdata( 2, 50, 1);%  model = ekozinec(data, struct('eps',0.01));%  figure; ppatterns(data); pline(model); %% See also %  PERCEPTRON, MPERCEPTRON, LINCLASS.%% About: Statistical Pattern Recognition Toolbox% (C) 1999-2003, Written by Vojtech Franc and Vaclav Hlavac% <a href="http://www.cvut.cz">Czech Technical University Prague</a>% <a href="http://www.feld.cvut.cz">Faculty of Electrical Engineering</a>% <a href="http://cmp.felk.cvut.cz">Center for Machine Perception</a>% Modifications:% 19-may-2004, VF% 3-may-2004, VF% 17-Sep-2003, VF% 17-Feb-2003, VF% 16-Feb-2003, VF% 21-apr-2001, V.Franc, created% get data dimensions[dim,num_data] = size(data.X);inx1=find(data.y==1);inx2=find(data.y==2);% Process input arguments% --------------------------if nargin < 2,  options = []; else options = c2s(options); endif ~isfield(options,'tmax'), options.tmax = inf; endif ~isfield(options,'eps'), options.eps = 0.01; endif ~isfield(options,'verb'), options.verb = 0; endif nargin < 3,  % creates init model  % --------------------------  inx1=find(data.y==1);   inx2=find(data.y==2);  model.W1 = data.X(:,inx1(1));  model.W2 = data.X(:,inx2(1));else  % take init model from input  %--------------------------------  model = init_model;endmodel.t = 0; model.exitflag = 0;% main loop%-----------------------------while model.exitflag == 0 & options.tmax > model.t,     model.t = model.t + 1;  dW = (model.W1 - model.W2);  norm_dW = norm( dW );    projx = data.X'*dW;    projx(inx1) = (projx(inx1) - model.W2'*dW)/norm_dW;  projx(inx2) = (-projx(inx2) + model.W1'*dW)/norm_dW;    [min_proj, min_inx] = min(projx);  % bound for separating or eps-optimal separating hyperplane   if options.eps < 0, bound = norm_dW/2; else bound=norm_dW-options.eps/2; end    if min_proj <= bound,    xt=data.X(:,min_inx);    % Updata - Kozinec's rule    if data.y(min_inx) == 1,      W1x = model.W1-xt;      k = min(1, dW'*W1x/ (W1x'*W1x));       model.W1 = model.W1 * (1 - k) + xt * k;    else      W2x = model.W2-xt;      k = min(1, -dW'*W2x / (W2x'*W2x));      model.W2 = model.W2 * (1 - k) + xt * k;    end    model.exitflag = 0;  else    model.exitflag = 1;  end  % print info  if options.verb == 1 & mod(model.t,100) == 0,     fprintf('iter %d: upper_bound = %f, lower_bound = %f, dif = %f\n', ...       model.t, norm_dW/2, min_proj/2, (norm_dW-min_proj)/2 );    end  endmodel.b = 0.5 * (model.W2'*model.W2 - model.W1'*model.W1);model.W = model.W1 - model.W2;model.margin = min_proj/2;model.fun = 'linclass';return;