disp('Running MATLAB script file DWDToyEg1.m') ;
%
%    Toy example, illustrating basics of 
%    Distance Weighted Discrimination and
%    Linear Support Vector Machine



%vipart = 4 ;
vipart = [0, 1, 2, 3, 4]  ;
                %          vector of indices of parts:
                %  0  - Raw Data
                %  1  - Bad plane
                %  2  - Better plane
                %  3  - Optimized DWD
                %  4  - Optimized SVM




for ipart = vipart ;

  fighand = figure(1) ;
  clf ;


  %  generate raw data
  %
  seed = 20943343 ;
  seed = 30658753 ;
  seed = 70937425 ;
  randn('seed',seed) ;
  n = 15 ;

  data1 = randn(1,n) ;
  data1 = [1.5 + data1 ; (0 + data1 + 0.6 * randn(1,n))] ;
  data2 = randn(1,n) ;
  data2 = [-1.5 + data2 ; (0 + data2 + 0.6 * randn(1,n))] ;

  bottom = min(min([data1 data2])) ;
  top = max(max([data1 data2])) ;
  range = top - bottom ;
  bottom = bottom - 0.02 * range ;
  top = top + 0.02 * range ;
  left = bottom - 0.15 * range ;
  right = top + 0.15 * range ;


  %  plot raw data
  %
  plot(data1(1,:),data1(2,:),'r+', ...
       data2(1,:),data2(2,:),'bo', ...
           'MarkerSize',8,'LineWidth',2) ;
    axis equal ;
    axis([left,right,bottom,top]) ;
    axis off ;
    title('Toy Discrimination Data','FontSize',24) ;
%    set(gca,'XTick',[]) ;
%    set(gca,'YTick',[]) ;

    text(left + 0.9 * (right - left), ...
         bottom + 0.3 * (top - bottom), ...
         'Class 1 (pos.)','Color','r','FontSize',18) ;

    text(left - 0.1 * (right - left), ...
         bottom + 0.7 * (top - bottom), ...
         'Class 2 (neg.)','Color','b','FontSize',18) ;

  pstr = 'DWDToyEg1' ;



  iplot = 0 ;
  if  ipart == 1  ;    %  do bad plane

    beta = 0 ;
    w = [0.1; -0.9] ;
    w = w / norm(w) ;

    iplot = 1 ;
    titstr = 'Toy Disc. Data, Poor Separating Plane' ;

      pstr = [pstr 'bad'] ;


  elseif  ipart == 2  ;    %  do good plane

    beta = 0 ;
    w = [0.45; -0.55] ;
    w = w / norm(w) ;

    iplot = 1 ;
    titstr = 'Toy Disc. Data, Better Separating Plane' ;

      pstr = [pstr 'good'] ;


  elseif  ipart == 3  ;    %  do optimized DWD
                           %  use lines modified from hdd1SM.m

    threshfact = 100 ;

    addpath '\Research\MatlabSoftware\HDD' -end ;   
    addpath '\Research\MatlabSoftware\HDD\SDPT3\Solver' -end ;   
    addpath '\Research\MatlabSoftware\HDD\SDPT3\Solver\mexexec' -end ;

    %  Compute median of pairwise distances squared between classes
    %
    vpwdist2 = [] ;
    for i = 1:n ;
      pwdist2 = sum((vec2matSM(data1(:,i),n) - data2).^2,1) ;
      vpwdist2 = [vpwdist2 pwdist2] ;
    end ;
    medianpwdist2 = median(vpwdist2) ;

    penalty = threshfact / medianpwdist2 ;

    [w,beta,residp,residn,totalviolation,X,y,Z] = sepold(data1,data2,penalty) ;

    [dr, dirvec] = hdd1SM([0; 0],data1,data2,11) ;
        %  use 0 as data, because don't care about result
        %  imeth = 11:    DWD (default param)

    iplot = 1 ;
    titstr = 'Toy Disc. Data, Optimized DWD' ;
    pstr = [pstr 'DWD'] ;


  elseif  ipart == 4  ;    %  do optimized SVM
                           %  use lines modified from hdd1SM.m

    threshfact = 1000 ;

    addpath \Research\MatlabSoftware\SVM\SVM1 -end ;

    d = size(data1,1) ;

    %  Compute median of pairwise distances squared between classes
    %
    vpwdist2 = [] ;
    for i = 1:n ;
      pwdist2 = sum((vec2matSM(data1(:,i),n) - data2).^2,1) ;
      vpwdist2 = [vpwdist2 pwdist2] ;
    end ;
    medianpwdist2 = median(vpwdist2) ;

    C = threshfact / medianpwdist2 ;
        %  threshfact "makes this large", 
        %  and 1 / medianpwdist2 "puts on correct scale"

    mtrainall = [data1 data2]' ;
    vclassall = [ones(1,n) -ones(1,n)]' ;
        %  since SVM wants data vectors as rows

    %  Run SVM algorithm
    %
    %  Using lines from program:  
    %      \Research\ComplexPopn\Old\PolyEmbed\PEegToyDisc1.m
    %

      ker = 'linear' ;
    [nsv,alpha,bias] = SVC(mtrainall,vclassall,ker,C) ;
        %  classification version of SVM

    %  Find normal vector to separating hyperplane
    %
    w = sum(vec2matSM(alpha,d) .* ...
                    mtrainall .* ...
            vec2matSM(vclassall,d),1) ;
        %  sum_i alpha_i x_i y_i

    %  Find intercept value of separating hyperplane
    %
    [temp,r] = max(alpha .* (vclassall == 1)) ;
        %  r is index of biggest alpha from class 1
    [temp,s] = max(alpha .*(vclassall == -1 )) ;
        %  s is index of biggest alpha from class -1
    beta = 0.5 * w * (mtrainall(r,:) + mtrainall(s,:))' ;

    w = w' / norm(w) ;
        %  make this a unit vector, and a column vector


    iplot = 2 ;
    titstr = 'Toy Disc. Data, Optimized SVM' ;
    pstr = [pstr 'SVM'] ;


  end ;    %  of ipart if-block



  if  iplot == 1  | ...
      iplot == 2  ;    %  then add normal vector, and separating hyperplane

    vcenter = beta * w ;
    hpdir = [w(2); -w(1)] ;
        %  vector in direction of hyperplane



    hold on ;


      title(titstr,'FontSize',24) ;

      plot([vcenter(1); vcenter(1) + 3 * w(1)], ...
           [vcenter(2); vcenter(2) + 3 * w(2)], ...
               'm-','LineWidth',6) ;
      text(vcenter(1) + 3 * w(1) + 0.3, ...
           vcenter(2) + 3 * w(2), ...
           'Normal Vector','Color','m','FontSize',18) ;

      plot([(vcenter(1) - 6 * hpdir(1)); (vcenter(1) + 6 * hpdir(1))], ...
           [(vcenter(2) - 6 * hpdir(2)); (vcenter(2) + 6 * hpdir(2))], ...
               'm--','LineWidth',6) ;
      text(vcenter(1) + 3 * hpdir(1) + 0.3, ...
           vcenter(2) + 3 * hpdir(2), ...
           'Separating Hyperplane','Color','m','FontSize',18) ;

      vproj1nv = data1' * w ;
      vproj2nv = data2' * w ;
            %  column vector of projections of data onto normal vector


      for i = 1:n ;

        vproj1hp = vcenter + (data1(:,i)' * hpdir) * hpdir ;
        vproj2hp = vcenter + (data2(:,i)' * hpdir) * hpdir ;
            %  projection of data point onto separating hyperplane

        if vproj1nv(i) >= beta ;    %  then data1 is on correct side
          plot([data1(1,i); vproj1hp(1,1)], ...
               [data1(2,i); vproj1hp(2,1)],'m-') ;
        else ;    %  then data1 is on wrong side
          plot([data1(1,i); vproj1hp(1,1)], ...
               [data1(2,i); vproj1hp(2,1)],'g--','LineWidth',3) ;
        end ;

        if vproj2nv(i) <= beta ;    %  then data2 is on correct side
          plot([data2(1,i); vproj2hp(1,1)], ...
               [data2(2,i); vproj2hp(2,1)],'m-') ;
        else ;    %  then data1 is on wrong side
          plot([data2(1,i); vproj2hp(1,1)], ...
               [data2(2,i); vproj2hp(2,1)],'g--','LineWidth',3) ;
        end ;

        text(left + 0.0 * (right - left), ...
             bottom + 0.93 * (top - bottom), ...
             'Residuals, r_i','Color','m','FontSize',15) ;

        text(left + 0.0 * (right - left), ...
             bottom + 0.86 * (top - bottom), ...
             'Violations, v_i','Color','g','FontSize',15) ;

      end ;


      if  iplot == 2  ;    %  then add additional SVM graphics

        marg = min(abs([vproj1nv; vproj2nv])) 
        flag = abs(marg - abs([vproj1nv; vproj2nv])) < 10^(-2) ;
        sum(flag)

        plot([(vcenter(1) + marg * w(1) - 6 * hpdir(1)); ...
                  (vcenter(1) + marg * w(1) + 6 * hpdir(1))], ...
             [(vcenter(2) + marg * w(2) - 6 * hpdir(2)); ...
                  (vcenter(2) + marg * w(2) + 6 * hpdir(2))], ...
                 'k--') ;
        plot([(vcenter(1) - marg * w(1) - 6 * hpdir(1)); ...
                  (vcenter(1) - marg * w(1) + 6 * hpdir(1))], ...
             [(vcenter(2) - marg * w(2) - 6 * hpdir(2)); ...
                  (vcenter(2) - marg * w(2) + 6 * hpdir(2))], ...
                 'k--') ;

        plot(mtrainall(flag,1),mtrainall(flag,2), ...
                 'ks','MarkerSize',18,'LineWidth',3) ;

        text(left + 0.0 * (right - left), ...
             bottom + 0.79 * (top - bottom), ...
             'Support Vectors','Color','k','FontSize',15) ;


      end ;    %  of 2nd ipart if-block


    hold off ;

  end ;    %  of 1st ipart if-block





    orient landscape ;
  print('-dpsc', [pstr '.ps']) ;




end ;    %  of ipart loop






