%
% SVC: Construct a support vector classifier.
%
% Constructs a support vector classifier based on the data points 'xo'
% (stored rowwise) and labels 'yo', which should be either -1 or 1.
%
% Returns:
%   'a'  a vector of weights for each point in 'xo'
%   'wo'
%   'bo' parameters in f(z) = wo.z + bo
%   'xs' a vector of support vectors
%   'ys' a vector of support vector labels
%
% NOTE:
%
% - Various options are not parameters; they can only be adjusted in this
%   file. See the list 'Parameters' below.
%
% - You'll need either the Matlab minimization package, which contains 'qp.m'
%   or some other quadratic minimization problem solver. Matlab's 'qp.m'
%   is not very stable, so it might give problems.
%
% - This is a rough, non-optimized version. Its only goal is to help
%   understand the algorithm.
%
% (C) 1996 David Tax & Dick de Ridder
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands
% Dick de Ridder, D.deRidder@ph.tn.tudelft.nl

function [ a, wo, bo, xs, ys ] = svc (xo, yo, options)

  % Some constants.

  false          = 0;
  true           = 1;

  % Methods.

  linear        = 0;
  polynomial    = 1;
  potential      = 2;
  radial_basis  = 3;
  neural        = 4;

  % Parameters.

  num_points        = size (xo, 1);
  input_dimension    = size (xo, 2);

  a_crit             = 1.0e-8;      % The criterion for considering a == 0.
  upper_bound        = 1.0e+20;    % Upper bound for a(i)'s.
  batch_size        = num_points;  % The batch size.
  max_retries        = 10;          % Maximum number of retries before the

  if nargin < 3

    % algorithm is restarted.
    feature_dimension  = 5;

    % The feature space dimensionality.
    method            = polynomial;     % Dot product convolution function.

    C  = 1.0e+2;        % Weighing factor of errors made.

    % Radial basis

    rb_sigma          = 1.0;

    % Potential

    po_sigma          = 1.0;

    % Neural

    nn_v              = 1.0e-1;
    nn_c              = 1.0;
    dographics        = true;
    dographics3d      = true;
  else
    method            = options(1);     % Dot product convolution function.
    feature_dimension = options(2);
    C                 = options(3);
    rb_sigma          = options(4);
    po_sigma          = options(5);
    nn_v              = options(6);
    nn_c              = options(7);
    dographics        = options(8);
    dographics3d      = options(9);
  end

  % Flags.

  dodebug      = false;
  throwaway    = false;

  % QP

  negdef      = 0;
  normalize    = 1;

  % Normalize: subtract mean.
  mx = mean (xo);
  for i = 1:input_dimension
    xo(:, i) = xo(:, i) - mx(i);
  end

  % Graphics: define the axes.

  mx = max (max (abs (xo))) * 1.05;
  ax = [-mx mx -mx mx];

  % We have four datasets:
  %    xorg   - the original dataset.
  %    xo     - the 'current' dataset, with processed points thrown
  %             out (keeps shrinking).
  %    x      - the current batch being processed: support vectors
  %             remaining from a previous cycle + some points randomly
  %             chosen from xo.
  %    xs     - the set of current support vectors.

  xorg = xo;
  yorg = yo;

  support_size   = 0;  xs = []; ys = []; a = [];

  org_num_points = num_points;

  done            = false;
  retries         = 0;

  % As long as there are unprocessed points....

  while ((size (xo, 1) > 0) & (done == false))

    % Prepare a new batch by copying the previously found
    % support vectors + a set of points randomly chosen from
    % xo (the set of remaining unprocessed points) into
    % x and y.

    points_left = size (xo, 1);

    num_points = min (batch_size + support_size, ...
                      points_left + support_size);
    r = randperm (points_left);
    r = r (1, 1:min (batch_size, points_left));

    x = [xs(1:support_size, :); xo(r, :)];
    y = [ys(1:support_size)   ; yo(r)   ];

    % Initialize D.

    D = zeros (num_points, num_points);
    for i = 1:num_points
      for j = i:num_points
        if (method == linear)
          D (i, j) = y(i) * y(j) * (x(i,:) * x(j,:)');
        elseif (method == polynomial)
          D (i, j) = y(i) * y(j) * ((x(i,:) * x(j,:)' + 1)^feature_dimension);
        elseif (method == potential)
          D (i, j) = y(i) * y(j) * exp (- vectdist (x(i,:), x(j,:)') / po_sigma);
        elseif (method == radial_basis)
          D (i, j) = y(i) * y(j) * exp (- (vectdist (x(i,:), x(j,:)') / rb_sigma)^2);
        elseif (method == neural)
          D (i, j) = y(i) * y(j) * tanh (nn_v * x(i,:) * x(j,:)' + nn_c);
        else
          disp ('Error: unknown dot product convolution.')
          return
        end
        D (j, i) = D (i, j);
      end
    end

    % Now check whether D is positive semi-definite, as it is supposed
    % to be. Thanks go to A.J. Quist of TWI-SSOR for the pd_check routine.

    if (pd_check (D) ~= 1)
      i = -50;
      while (pd_check (D + (10.0^i) * eye (num_points)) ~= 1)
        i = i + 1;
      end

      fprintf ('Adding I * 10^%d to D\n', i)

      D = D + (10.0^(i)) * eye(num_points);
    end

    D = [D                    zeros(num_points,1)  ; ....
         zeros(1,num_points)   1/C];

    % Initialize a(i) randomly

    rand ('seed', sum (100 * clock));
    a = [ 0.1 * rand(1, num_points) 1];

    % Initialize graphics and freeze the axes.

    if (dographics)

      % Init.

       figure(1);
       clf;
      axis (ax); axis ('square');

      hold on;

      % Plot the points

      plot (xorg (find (yorg < 0), 1), xorg (find (yorg < 0), 2), 'y+');
      plot (xorg (find (yorg > 0), 1), xorg (find (yorg > 0), 2), 'y*');
      plot (xo   (find (yo   < 0), 1), xo     (find (yo   < 0), 2), 'g+');
      plot (xo   (find (yo   > 0), 1), xo     (find (yo   > 0), 2), 'g*');
      plot (x    (find (y    < 0), 1), x    (find (y    < 0), 2), 'co');
      plot (x    (find (y    > 0), 1), x    (find (y    > 0), 2), 'co');

      if (dodebug)
        for i = 1:num_points
          text (x(i, 1) + 0.1, x(i, 2) + 0.1, num2str (i));
        end
      end

      drawnow

      if (dodebug)
        wait ('Press any key');
      end

    end

    % Minimization procedure initialization:
    % 'qp' minimizes:   0.5 x' D x + f' x
    % subject to:       Ax <= b
    %
    % Here: one constraint is coded into the lower
    % bound of the solution (i.e., a(i) >= 0, eqn. 16),
    % the other is given in f:
    % Ax <= b with A = vector containing y(i)'s,
    % x = the solution to be found (a(i)'s) and
    % b = 0 gives a'y = 0 (eqn. 17). So we set f to
    % (-1 -1 ... -1)'.
    %
    % NB: Since we want to *maximize* eqn. 26,
    %     we minimize its negative - that's why
    %     f is negative.

    f = [-ones(num_points, 1); 0];

    % The first row in A is the equality constraint
    % ay = 0. The second to last rows are the constraints
    % given in eqn. 29: all a(i) should be smaller than
    % delta....

    A = [y'              0                     ; ...
         eye(num_points) -ones(num_points, 1)  ];
    b = zeros (num_points + 1, 1);

    % One of our constraints: each a(i) should be >= 0 (the
    % lower bound, lb). The upper bound is an arbitrary
    % number.

    lb   = zeros (num_points + 1, 1);
     ub   = [ upper_bound * ones(num_points + 1, 1)];

    % Initial guess. This doesn't seem to have an
    % enormous effect.

    rand ('seed', sum( 100 * clock));
    p = 0.5 * rand (num_points + 1, 1);

    % Call the minimization routine. Note: this routine
    %  has some built-in safety checks and returns []
    % when it decides something is definitely wrong.
    %
    % The last 1 in the call means that only the first
    % constraint is an equality constraint. See 'qp'.

    if (exist ('qld') == 3)
      a = qld (D, f, -A, b, lb, ub, p, 1);
    else
      a = qp (D, f, A, b, lb, ub, p, 1, 1, negdef, normalize)';
    end

    if (~isempty(a))
      delta = a(num_points + 1);
      a(num_points + 1) = [];
    else
      disp ('Geen convergentie')
    end

    % If this minimization went OK, find the
    % support vectors and throw all points from this batch
    % out of the pool.
    %
    % OK means: There is at least one support vector
    % (i.e., a(i) > a_crit for some i) and we do not have
    % a vector of NaN's or an empty vector.

    if ((prod(a <= (a_crit * ones (size (a,1),1))) == 0) & ...
        (isempty(find(isnan(a)==1))) & ...
        (~isempty(a)) & ...
        (delta + 0.001 >= max(a)))

      % Throw the set of processed points out of the
      % set of remaining points.

      points_left = points_left - num_points;

      xo(r,:) = [];
      yo(r)   = [];

      % Find all support vectors in the current
      % batch. Plot 'real' support vectors in red,
      % overlapping vectors (i.e., errors) in magenta.

      support_size = 0;
      for i=1:num_points
        if (a(i) > a_crit)
          if (dographics)
            if (a(i) < delta-0.0001)
              plot (x(i,1), x(i,2), 'ro');
            else
              plot (x(i,1), x(i,2), 'mo');
            end
          end
          support_size = support_size + 1;
            xs(support_size, :) = x(i, :);
            ys(support_size, :) = y(i);
        end
      end

      % Optimal hyperplane: calculate w.

      if (method == linear)
        wo = zeros (feature_dimension, 1);
          for i = 1:input_dimension
            for j = 1:num_points
              wo(i) = wo(i) + a(j)*y(j)*x(j,i);
           end
          end
      end

      % Calculate b, the offset of the classifying function.
      % Note that this only works for 'real' support vectors,
      % i.e. non-errors. We believe that we can discriminate
      % between these on the basis that errors have a(i) == delta,
      % unlike real support vectors. This is a bit tricky, however.

      real_sv = find ((a > a_crit) & (a < delta - 0.0001));
      if (isempty(real_sv))
        real_sv = find (a > a_crit);
      end
      bo = zeros (size (real_sv));

      for j = 1:size(real_sv,1)
        r = real_sv(j);
        g = 0;
        for i = 1:num_points
          if (method == linear)
            g = g + y(i) * a(i) * (x(r,:) * x(i,:)');
          elseif (method == polynomial)
            g = g + y(i) * a(i) * ((x(r,:) * x(i,:)' + 1)^feature_dimension);
          elseif (method == potential)
            g = g + y(i) * a(i) * exp (- vectdist (x(r,:), x(i,:)') / po_sigma);
          elseif (method == radial_basis)
            g = g + y(i) * a(i) * exp (- (vectdist (x(r,:), x(i,:)') / rb_sigma)^2);
          elseif (method == neural)
            g = g + y(i) * a(i) * tanh (nn_v * x(r,:) * x(i,:)' + nn_c);
          else
            disp ('Error: unknown dot product convolution.')
            return
          end
        end
        bo(j) = y(r) - g;
        b = bo;
      end
      b=mean(bo); bo=b;

      % Draw the classifier and its margins.

      if (dographics & (method == linear))
        h=line ([-mx mx], [((-wo(1)*(-mx))-b  )/wo(2) ((-wo(1)*mx)-b)/wo(2)]);
        set(h,'Color','w');
         h=line ([-mx mx], [((-wo(1)*(-mx))-b-1)/wo(2) ((-wo(1)*mx)-b-1)/wo(2)]);
        set(h,'Color','r');
         h=line ([-mx mx], [((-wo(1)*(-mx))-b+1)/wo(2) ((-wo(1)*mx)-b+1)/wo(2)]);
        set(h,'Color','r');
      end

      % Throw out all points which lie beyond the margin. Quote the Raven:
      % 'These will become a support vector .... NeverMore.'
      %
      % This does not work if we do not consider points that have been
      % rejected. Solution: build a set of the points in the *entire*
      % dataset that do not satisfy constraint eqn. 10 (David).

      if (throwaway)
        r = zeros (size (xorg, 1), 1);
        for i = 1:size(r,1)
           g = 0;
           for j = 1:num_points
             if (method == linear)
               g = g + y(j) * a(j) * (xorg(i,:) * x(j,:)');
             elseif (method == polynomial)
               g = g + y(j) * a(j) * ((xorg(i,:) * x(j,:)' + 1)^feature_dimension);
             elseif (method == potential)
               g = g + y(j) * a(j) * exp (- vectdist (xorg(i,:), x(j,:)') / po_sigma);
             elseif (method == radial_basis)
               g = g + y(j) * a(j) * exp (- (vectdist (xorg(i,:), x(j,:)') / rb_sigma)^2);
            elseif (method == neural)
              g = g + y(j) * a(j) * tanh (nn_v * xorg(i,:) * x(j,:)' + nn_c);
             else
               disp ('Error: unknown dot product convolution.')
               return
             end
           end
          r(i) = yorg(i) * (g + b);
        end

        xo = xorg (find (r < 1),:);
        yo = yorg (find (r < 1),:);

      end

      %  !%#(^%&!^#$ We now get the possibility of having a
      %  support vector in set xs *and* in xo because r < 1.
      % (soft margin classifier). Therefore: kick any
      %  points in xs out of xo.

      k = size (xo,1);
      for j = 1:support_size
        i = 1;
        while i <= k
          if (vectdist(xo(i,:),xs(j,:)') == 0)
            xo(i,:) = [];
            yo(i,:) = [];
            k = k - 1;
          else
            i = i + 1;
          end
        end
      end

      % Clip a (throw away non-support vector values).

      a  = a(find (a > a_crit));

      % Clip the set of support vectors.

      xs = xs(1:support_size, :);
      ys = ys(1:support_size);

      if (dodebug)
         wait ('Press any key');
      end

      s = fprintf ('%d vectors to go, %d support \n', ...
                    size (xo, 1), support_size);

      retries = 0;

    else

      retries = retries + 1;

      % If all else fails, start over....

      if (retries > max_retries)
        if (dodebug)
          disp ('Re-starting algorithm.')
        end
        support_size = 0;  xs = []; ys = []; a = [];
        xo = xorg;
        yo = yorg;
        retries = 0;
      end

    end

    if (dodebug)
      disp ('Number of points left to process:')
      size(xo,1)
    end

  end

  if (~dodebug)
    clf;
    hold on;
  end

  % Plot the classifier, all datapoints and mark the support vectors.

  if (dographics)

    if (method == linear)
      h = line ([-mx mx], [((-wo(1)*(-mx))-b  )/wo(2) ((-wo(1)*mx)-b )/wo(2)]);
        set (h, 'Color', 'w');
       h = line ([-mx mx], [((-wo(1)*(-mx))-b-1)/wo(2) ((-wo(1)*mx)-b-1)/wo(2)]);
       set (h, 'Color', 'r');
       h = line ([-mx mx], [((-wo(1)*(-mx))-b+1)/wo(2) ((-wo(1)*mx)-b+1)/wo(2)]);
        set (h, 'Color', 'r');
    end

    plot (xorg(find(yorg<0),1), xorg(find(yorg<0),2), 'y+');
    plot (xorg(find(yorg>0),1), xorg(find(yorg>0),2), 'y*');
    plot (xs(find(a<delta-0.0001),1), xs(find(a<delta-0.0001),2), 'ro');
    plot (xs(find(a>=delta-0.0001),1), xs(find(a>=delta-0.0001),2), 'mo');

    axis (ax); axis ('square');
    drawnow;

  end

  if (dographics3d & (input_dimension == 2))

    fprintf ('Plotting function in 3D\n')

    figure (2);

    [fx,fy] = meshgrid (-mx:.03:(mx*1.01),-mx:.03:(mx*1.01));
    for i = 1:size(fx, 1)
      for j = 1:size(fx, 2)
        fz(i, j) = b;
        for k = 1:support_size
          if (method == linear)
            fz(i, j) = fz(i, j) + ys(k) * a(k) * ...
            ([fx(i,j) fy(i,j)] * xs(k,:)');
          elseif (method == polynomial)
            fz(i, j) = fz(i, j) + ys(k) * a(k) * ...
            (([fx(i,j) fy(i,j)] * xs(k,:)' + 1)^feature_dimension);
          elseif (method == potential)
            fz(i, j) = fz(i, j) + ys(k) * a(k) * ...
            exp (- vectdist ([fx(i,j) fy(i,j)], xs(k,:)') / po_sigma);
          elseif (method == radial_basis)
            fz(i, j) = fz(i, j) + ys(k) * a(k) * ...
            exp (- (vectdist ([fx(i,j) fy(i,j)], xs(k,:)') / rb_sigma)^2);
          elseif (method == neural)
            fz(i, j) = fz(i, j) + ys(k) * a(k) * ...
            tanh (nn_v * [fx(i,j) fy(i,j)] * xs(k,:)' + nn_c);
          else
            disp ('Error: unknown dot product convolution.')
            return
          end
        end
      end
    end

    figure (2);
    clf;
    hold on;
    surf (fx, fy, fz);
    shading flat;
    colormap('default'); colormap(colormap/2);
    contour3 (fx, fy, fz, [0 0], 'k');
    contour3 (fx, fy, fz, [-1 -1], 'w-.');
    contour3 (fx, fy, fz, [ 1  1], 'w-.');
    plot3 (xorg(find(yorg<0),1), xorg(find(yorg<0),2), ...
            10*max(max(fz))*ones(size(find(yorg<0),1),1), 'w+');
     plot3 (xorg(find(yorg>0),1), xorg(find(yorg>0),2), ...
            10*max(max(fz))*ones(size(find(yorg>0),1),1), 'w*');
    plot3 (xs(:,1), xs(:,2), 100*ones(support_size,1), 'ko');
    colorbar
    axis ([min(min(fx)) max(max(fx)) min(min(fy)) max(max(fy))]);
    axis ('square');
  end

return

