demo_linclass.m

很好的matlab模式识别工具箱

    tmp_model = ekozinec(sets,tmp_options, init_model);
   else
     tmp_model = ekozinec(sets,tmp_options);
   end
   
   if iter==-1,
     if tmp_model.last_update, alpha=sets.X(:,tmp_model.last_update); 
     else alpha=tmp_model.W; end
   else
     alpha = tmp_model.W;
   end
   solution=tmp_model.exitflag;
   tplus1 = tmp_model.t;
   theta = -tmp_model.b;

 case 3
%    [alpha,theta,solution,tplus1]=ekozinec(sets.X,sets.I,epsil,iter,t,...
%                                           alpha,theta);      
  if iter==-1, tmp_options.tmax = -1; else tmp_options.tmax=t+iter; end
   if t~=0, 
    init_model.t=t;
    init_model.W=alpha;
    init_model.b=-theta;
    tmp_options.eps=epsil;
    tmp_model = ekozinec(sets,tmp_options,init_model);
   else
    tmp_options.eps=epsil;
    tmp_model = ekozinec(sets,tmp_options);
   end
   
   if iter==-1,
     if tmp_model.last_update, alpha=sets.X(:,tmp_model.last_update); 
     else alpha=tmp_model.W; end
   else
     alpha = tmp_model.W;
   end
   solution=tmp_model.exitflag;
   tplus1 = tmp_model.t;
   theta = -tmp_model.b;

 case 4
    tmp_options.ker='linear';
    tmp_options.C = 1e12;
%    [alpha,theta,solution]=svmmot(sets.X,sets.I);
    tmp_model=smo(sets,tmp_options);
%    if tmp_model.exitflag <=0, solution=-1; else solution=1; end
    if tmp_model.trnerr > 0, solution = -1; else solution=1; end
    theta=-tmp_model.b;
    alpha=tmp_model.sv.X*tmp_model.Alpha(:);
    tplus1=1;
 case 5
%    [alpha,theta,solution,tplus1]=psum(sets.X,sets.I,10,iter,t,alpha,theta);
% case 6
%    [alpha,theta,solution,tplus1]=psumv(sets.X,sets.I,iter,t,alpha,theta);
% case 7
%    [alpha,theta,solution,tplus1]=navara1(sets.X,sets.I,iter,t,alpha,theta);
% case 8
%    [alpha,theta,solution,tplus1]=navarah1(sets.X,sets.I,iter,t,alpha,theta);
% case 9
%    [alpha,theta,solution]=simplex(sets.X,sets.I);
%    if solution==0, solution=-1; end
%    tplus1=1;

end

return


function []=setaxis(handle,rect)                                                
% function []=setaxis(handle,rect)                                              
%                                                                               
% SETAXIS sets scaling for the x- and y-axes                                    
%   on the plot with a given handle.                                            
%                                                                               
% See also AXIS.                                                                
%                                                                               
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac         
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz                
% Written Vojtech Franc (diploma thesis) 02.01.2000                             
% Modifications                                                                 
                                                                                
set(handle,'XLim',rect(1:2));                                                   
set(handle,'YLim',rect(3:4));                                                   
                                                                                
if size(rect,2)>=6,                                                             
   set(handle,'ZLim',rect(5:6));                                                
end                                                                             
                                                                                
return;     

%========================================                                       
function [rect]=getaxis(handle)                                                 
% function [rect]=getaxis(handle)                                               
%                                                                               
% GETAXIS returns a row vector containing the scaling for                       
%   the plot with a given handle.                                               
%                                                                               
% See also AXIS.                                                                
%                                                                               
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac         
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz                
% Written Vojtech Franc (diploma thesis) 02.01.2000                             
% Modifications                                                                 
% 24. 6.00 V. Hlavac, comments polished.                                        
                                                                                
rect=[get(handle,'XLim'),get(handle,'YLim'),get(handle,'ZLim')];                
                                                                                
return;      

function [win]=cmpwin(mins,maxs,xborder,yborder)                                
%                                                                               
%  [win]=cmpwin(mins,maxs,xborder,yborder)                                      
%                                                                               
% CMPWIN computes appropriate size of the axes.                                 
%                                                                               
                                                                                
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac         
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz                
% Written Vojtech Franc (diploma thesis) 02.01.2000                             
% Modifications                                                                 
% 24. 6.00 V. Hlavac, comments polished.                                        
                                                                                
dx=max( (maxs(1)-mins(1)), 1 )*xborder;                                         
dy=max( (maxs(2)-mins(2)), 1 )*yborder;                                         
                                                                                
%x1=floor(mins(1)-dx);                                                          
%x2=ceil(maxs(1)+dx);                                                           
%y1=floor(mins(2)-dy);                                                          
%y2=ceil(maxs(2)+dx);                                                           
x1=(mins(1)-dx);                                                                
x2=(maxs(1)+dx);                                                                
y1=(mins(2)-dy);                                                                
y2=(maxs(2)+dx);                                                                
                                                                                
win=[x1 x2 y1 y2];  

return;


function [x1,y1,x2,y2,in]=cliplin1(alpha,theta,window)                          
% [x1,y1,x2,y2,in]=cliplin1(alpha,theta,window)                                 
%                                                                               
% CLIPLIN1 clips the line given by the equation alpha*x=theta along             
%   the window. It returns two points on the border of the window.              
%   If the line is in the window then the argument is equal to 1                
%   else it returns 0.                                                          
%                                                                               
% See also CLIPLIN2.                                                            
%                                                                               
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac         
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz                
% Written Vojtech Franc (diploma thesis) 20.10.1999, 23.12.1999                 
% Modifications                                                                 
% 24. 6.00 V. Hlavac, comments polished.                                        
                                                                                
minx=window(1);                                                                 
maxx=window(2);                                                                 
miny=window(3);                                                                 
maxy=window(4);                                                                 
                                                                                
x=zeros(4,1);                                                                   
y=zeros(4,1);    
if alpha(1)==0,                                                                 
   if alpha(2)~=0,                                                              
      x1=minx;                                                                  
      y1=theta/alpha(2);                                                        
      x2=maxx;                                                                  
      y2=y1;                                                                    
      in=1;                                                                     
   else                                                                         
      % if alpha == 0 then it means the bad input.                              
      x1=0;                                                                     
      y1=0;                                                                     
      x2=0;                                                                     
      y2=0;                                                                     
      in=0;                                                                     
   end                                                                          
elseif alpha(2)==0,                                                             
   x1=theta/alpha(1);                                                           
   y1=miny;                                                                     
   x2=x1;                                                                       
   y2=maxy;                                                                     
   in=1;                                                                        
 else      
     y(1)=maxy;                                                                   
   x(1)=(theta-alpha(2)*y(1))/alpha(1);                                         
   y(2)=miny;                                                                   
   x(2)=(theta-alpha(2)*y(2))/alpha(1);                                         
                                                                                
   x(3)=maxx;                                                                   
   y(3)=(theta-alpha(1)*x(3))/alpha(2);                                         
   x(4)=minx;                                                                   
   y(4)=(theta-alpha(1)*x(4))/alpha(2);                                         
                                                                                
   j=0;                                                                         
   for i=1:4,                                                                   
      if x(i) <= maxx & x(i) >= minx & y(i) <= maxy & y(i) >= miny,             
         if j==0,                                                               
            j=j+1;                                                              
            x1=x(i);                                                            
            y1=y(i);                                                            
         elseif j==1,                                                           
            j=j+1;                                                              
            x2=x(i);                                                            
            y2=y(i);                                                            
         end                                                                    
      end                                                                       
   end                                                                          
                                                                                
   if j<2,                                                                      
      x1=0;                                                                     
      y1=0;                                                                     
      x2=0;                                                                     
      y2=0;                                                                     
      in=0;                                                                     
   else                                                                         
      in=1;                                                                     
   end                                                                          
end % elseif alpha(2)==0                                                        
     

function []=clrchild(handle)                                                    
% function []=clraxis(handle)                                                   
%                                                                               
% CLRCHILD clears children of an object with the given handle.                  
%                                                                               
% See also DELETE.                                                              
%                                                                               
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac         
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz                
% Written Vojtech Franc (diploma thesis) 02.01.2000                             
% Modifications                                                                 
% 24. 6.00 V. Hlavac, comments polished.                                        
                                                                                
delete(get(handle,'Children'));                                                 
                                                                                
return;    

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
%  [Koz73] with eps-optimality stopping condition ([SH10], Chapter 5).
%  It transforms the binary-class problem to the single-class problem.
%  The nearest point from convex hull of transformed data which
%  determines the optimal separating hyperplane is found by the 
%  Kozinec's algorithm.
%  
%  model = ekozinec(data) finds eps-optimal separating hyperplane for
%   linearly separable data.
%
%  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 Euclidien norm of the
%        normal vector (default 0.01).
%    .tmax [1x1]... maximal number of iterations.
%  
%   If tmax==-1 then it only returns index (model.last_update)
%   of data vector which should be used by the algorithm for updating
%   the linear rule in the next iteration.
%  
%  model = ekozinec(data,options,init_model) specifies initial model
%   which must contain:
%    .W [dim x 1] ... normal vector.
%    .b [1x1] ... bias of hyperplane.
%    .t [1x1] (optional) ... iteration number.
%
% 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 [1x1] Precision of found hyperplane (see above).
%   .tmax [1x1] Maximal number of iterations (default tmax=inf).
%     If tmax==-1 then it does not perform any iteration but returns only 
%     index of the point which should be used to update linear rule.
%  
%  init_model [struct] Initial model; must contain items
%    .W, .b and .t (see below).
%
% Output:
%  model [struct] Found linear classifier:
%   .W [dim x 1] Normal vector of hyperplane.
%   .b [1x1] Bias of hyperplane.
%  
%   .exitflag [1x1] 1 ... eps-optimality condition satisfied
%                   0 ... number of iterations exceeded tmax.
%   .t [int] Number of iterations.
%   .last_update [d x 1] Index of the last point used for update.
%
% Example:
%  data = genlsdata( 2, 50, 1);
%  model = ekozinec(data)
%  figure; ppatterns(data); pline(model); 
%
% See also 
%

% 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:
% 3-may-2004, VF
% 17-sep-2003, VF
% 16-Feb-2003, VF
% 21-apr-2001, V.Franc, created

% get data dimensions
[dim,num_data] = size(data.X);

% Process input arguments 
% ------------------------------------
if nargin < 2, options = []; else options = c2s( options ); end

if ~isfield(options,'tmax'), options.tmax = inf; end
if ~isfield(options,'eps'), options.eps = 0.01; end

if nargin < 3, 
  % make inital model
  model.b = 0;
  if data.y(1)==1, W = [data.X(:,1);1]; else W = -[data.X(:,1);1]; end
else
  % get inital model from input
  model = init_model;
  W = [model.W; model.b];
end
if ~isfield( model,'t'), model.t = 0; end

model.exitflag = 0;
model.last_update = 0;

% Add one constant coordinates to the data and swap                             
% points from the second class along the origin.
% ---------------------------------------------------- 
data.X = [data.X; ones(1,num_data)];
dim=dim+1;
inx = find(data.y==2);
data.X(:,inx) = -data.X(:,inx);

if options.tmax == -1,
  % return index of point which should be used to update linear rule
  %----------------------------------------------------------------------
  norm_W = norm(W);
  [min_proj,min_inx]=min( data.X'*W/norm_W );

  % bound for separating or eps-optimal separating hyperplane 
  if options.eps < 0, bound = norm_W/2; else bound = norm_W - options.eps; end

  model.last_update = min_inx;
  if min_proj <= bound, model.exitflag = 0; else model.exitflag = 1; end
else 

  % main loop
  %----------------
  while model.exitflag == 0 & options.tmax > model.t,

    model.t = model.t + 1;

    norm_W = norm(W);
    [min_proj,min_inx]=min( data.X'*W/norm_W );

    % bound for separating or eps-optimal separating hyperplane 
    if options.eps < 0, bound = 0; else bound=norm_W - options.eps; end
    
    if min_proj <= bound,
      model.last_update = min_inx;
      xt=data.X(:,min_inx);
    
      k=min(1,(W-xt)'*W/((W-xt)'*(W-xt)) );
      W = W*(1-k) + xt*k;
      
      model.exitflag = 0;
    else
      model.exitflag = 1;
    end
  end
  
end

model.W = W(1:dim-1);
model.b = W(dim);
model.fun = 'linclass';

return;