function [alphas,solution,t]=fisherk(X,J,K,tmax,t,alphas)
% FISHERK learns the Fisher classifier using Kozinec's rule.
%  [alphas,solution,t]=fisherk(X,J,K,tmax,t,alphas)
%
% FISHERK This algorithm finds the Fisher`s classifier
%   by the use of the modified Kozinec`s algorithm.
%   This task is equivalent to finding solution of
%   the array of linear nonequalities. The goal is to
%   find such alphas(:,i), i=1,2,...K that hold
%       alphas(:,i)'*X(:,k) > alphas(:,j)'*X(:,k)  for J(k)=i, i ~= j
%
%   This algorithm works iteratively and find solution in finite
%   number of steps if the solution exists. 
%
% Input:
%  X [NxM] matrix containing M training points in N-dimensional
%     feture space. X=[x1,x2,...xM].
%  J [1xM] vector containing M integer class labels for
%     each point from X. Possible are integer values from 1 to K.
%  K [1x1] is number of classes.
%  tmax [1x1] is upper limit of number of algorithm steps.
%
%  t [1x1], alphas [NxK] if these arguments enter function then
%     the algorithm starts up from the state they define.
%     The argument t is an initial step number and matrix 'alphas' 
%     contains an initial solution.
%
% Output:
%  alphas [NxK] contains found solution. Vector alphas(:,i) coresponds 
%     the to i-th class.
%  solution [1x1] is equal to 1 if solution was found.
%                 is equal to 0 if was not found.
%  t [1x1] number of iterations.
%
% See also FISHERP, FISHDEMO, PERCEPTR, KOZINEC.
%

% 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
% 26-June-2001, V.Franc, comments improved.
% 24. 6.00 V. Hlavac, comments polished.

% get dim. and # of points
N=size(X,1);
NX=size(X,2);

% default setting
if nargin < 5,
   t=0;
end
if nargin < 4,
   tmax=inf;
end

% STEP (1)
if t==0,
   alphas=zeros(N,K);
   for i=1:K,
      index=find(J==i);
      if isempty(index)~=1,
         % get the first found vector
         alphas(:,i)=X(:,index(1));
      end
   end
end

% iterate until solution is found
solution = 0;
while solution == 0 & tmax > 0,
   tmax = tmax-1;
   solution=1;

   for i=1:NX,
       b=alphas'*X(:,i);
      [b,k]=max(b);
      j=J(i);

      if k ~= j,
         % adjust alpha

         % transformation
         alpha=[];
         x=zeros(1,K*N);
         for m=1:K,
            alpha=[alpha,alphas(:,m)'];
            if m==k,
               x((m-1)*N+1:m*N)=-X(:,i)';
            end
            if m==j,
               x((m-1)*N+1:m*N)=X(:,i)';
            end
         end

         % l = min(1, argmin | alpha*(1-l) + x*l |)
         l=-min(1,((x-alpha)*alpha')/((x-alpha)*(x-alpha)'));

         alphas(:,j)=alphas(:,j)+X(:,i)*l/(1-l);
         alphas(:,k)=alphas(:,k)-X(:,i)*l/(1-l);
         alphas=alphas*(1-l);
         t=t+1;
         solution=0;
         break;
      end
   end
end