disp('Running MATLAB script file DWD1figB.m') ;
%
%    Makes Figure B for the Distance Weighted Discrimination paper
%
%    Copied from SVMeg3.m, which:
%
%    Investigates  50-d example of use of the
%    Linear Support Vector Machine
%
%    Major change is to use mu0 = +- 2.2,
%    as in:    hdd1sim1.m


iout = 2 ;      %  1 - Movie of increasing dimension
                %  2 - Figure for paper



imeth = 5 ;
                 %  5 - FLD, but common mean covariance est.




moviefps = 5 ;
moviecstr = 'Cinepak' ;
nrepeat = 2 ;

dsave = 39 ;
    %  dimension to save for paper


left = -4 ;
right= 4 ;
bottom = -2 ;
top = 2 ;

n1 = 20 ;
n2 = 20 ;



seed = 7543870734 ;
randn('state',seed) ;


mu0 = 2.2 ;
sig0 = 1 ;
mtraindat1 = mu0 + sig0 * randn(n1,1) ;
mtraindat2 = -mu0 + sig0 * randn(n2,1) ;

sig = 1 ;

mtraindat = [mtraindat1; mtraindat2] ;


vtrainclass = [ones(n1,1); -ones(n2,1)] ;
        %  class labels


h1 = bwsjpiSM(mtraindat1,[left; right]) ;
h2 = bwsjpiSM(mtraindat2,[left; right]) ;
h = exp(mean(log([h1; h2]))) ;
  paramstruct = struct('vxgrid',[left; right], ...
                       'vh',h) ;
[d1kde1,kdexgrid] = kdeSM(mtraindat1,paramstruct) ;
[d2kde1,kdexgrid] = kdeSM(mtraindat2,paramstruct) ;
top1 = max([d1kde1; d2kde1]) ;
bottom1 = -0.05 * top1 ;
top1 = 1.05 * top1 ;
rhts1 = rand(n1,1) ;
rhts2 = rand(n2,1) ;
d1rhts1 = bottom1 + (.7 + .1 * rhts1) * (top1 - bottom1) ;
d2rhts1 = bottom1 + (.5 + .1 * rhts2) * (top1 - bottom1) ;
            %  random heights




fighand = figure(1) ;
clf ;
clear moviestruct ;



if imeth == 1 ;
  disctitstr = 'Fisher Linear Discrimination' ;
  outmethstr = 'FLD' ;
elseif imeth == 2 ;
  disctitstr = 'Linear SVM, C = 1000' ;
  outmethstr = 'SVM3' ;
elseif imeth == 3 ;
  disctitstr = 'Linear SVM, C = 10^{12}' ;
  outmethstr = 'SVM12' ;
elseif imeth == 4 ;
  disctitstr = 'Linear SVM, C = 10^{-6}' ;
  outmethstr = 'SVMm6' ;
elseif imeth == 5 ;
  disctitstr = 'FLD, Global Cov. based' ;
  outmethstr = 'FLDgc' ;
elseif imeth == 6 ;
  disctitstr = 'Mean Difference' ;
  outmethstr = 'MD' ;
end ;



%  Make dimension 1 plot
%
if iout == 1 ;
  ms = 10 ;
  lw = 3 ;
  lwkde = 2 ;
elseif iout == 2 ;
  ms = 6 ;
  lw = 1 ;
  lwkde = 1 ;
end ;
subplot(2,1,1) ;
plot(mtraindat1,zeros(n1,1),'r+',...
     mtraindat2,zeros(n2,1),'bo',...
         'MarkerSize',ms,'LineWidth',lw) ;
t1h = title([disctitstr ', dimension = 1'],'FontSize',15) ;
axis([left,right,bottom,top]) ;
hold on ;
  plot([left; right],[0; 0],'m--','LineWidth',2) ;
hold off ;
correctdir = 1 ;


subplot(2,2,3) ;
plot(kdexgrid,d1kde1,'r-', ...
     kdexgrid,d2kde1,'b-', ...
     [left; right],[0; 0],'m--','LineWidth',lwkde) ;
     axis([left,right,bottom1,top1]) ;
  t2h = title('True Projection','FontSize',15) ;
  hold on ;
    plot(mtraindat1,d1rhts1,'r+','MarkerSize',ms,'LineWidth',lw) ;
    plot(mtraindat2,d2rhts1,'bo','MarkerSize',ms,'LineWidth',lw) ;
  hold off ;


subplot(2,2,4) ;
plot(kdexgrid,d1kde1,'r-', ...
     kdexgrid,d2kde1,'b-', ...
     [left; right],[0; 0],'m-','LineWidth',lwkde) ;
     axis([left,right,bottom1,top1]) ;
  t3h = title('Est''d Projection','FontSize',15) ;
  hold on ;
    plot(mtraindat1,d1rhts1,'r+','MarkerSize',ms,'LineWidth',lw) ;
    plot(mtraindat2,d2rhts1,'bo','MarkerSize',ms,'LineWidth',lw) ;
  hold off ;



if iout == 1 ;    %  Increasing dimension movie

  moviestruct(1) = getframe(fighand) ;

  maxd = 50 ;

elseif iout == 2 ;

  maxd = dsave ;

end ;




for d = 2:maxd ;

  mtraindat1 = [mtraindat1, (sig * randn(n1,1))] ;
  mtraindat2 = [mtraindat2, (sig * randn(n2,1))] ;
  mtraindat = [mtraindat1; mtraindat2] ;
        %  Tack on a random column


  %  Run appropriate discriminator
  %
  if  imeth == 1  | ...
      imeth == 5  ;    %  then do some type of FLD

    vmean1fld = mean(mtraindat1',2) ;
    vmean2fld = mean(mtraindat2',2) ;
    vmeanfld = mean([vmean1fld vmean2fld],2) ;

    if imeth == 1 ;    %  then do proper FLD
      mresid1 = mtraindat1' - vec2matSM(vmean1fld,n1) ;
      mresid2 = mtraindat2' - vec2matSM(vmean2fld,n2) ;
    elseif imeth == 5 ;    %  then wnat global covariance
      mresid1 = mtraindat1' - vec2matSM(vmeanfld,n1) ;
      mresid2 = mtraindat2' - vec2matSM(vmeanfld,n2) ;
    end ;
    mresid = [mresid1 mresid2] ;

    mcov = cov(mresid') ;
        %  Get covariance matrix, transpose, since want 
        %               "coordinates as variables"
    mcovinv = pinv(mcov) ;
        %  pseudo-inverse

    fldv = mcovinv * (vmean1fld - vmean2fld) ;
        % Fisher Linear Discriminant Vector
    fldv = fldv / norm(fldv) ;
        %  make it a unit vector
    w = fldv' ;
        %  normal vector to separating plane


  elseif imeth == 6 ;    %  then do simple mean difference

    vmean1 = mean(mtraindat1',2) ;
    vmean2 = mean(mtraindat2',2) ;

    w = (vmean1 - vmean2) ;
    w = w / norm(w) ;
    w = w' ;
        %  normal vector to separating plane


  else ;    % then do an SVM

    if imeth == 2 ;
      C = 1000 ; 
    elseif imeth == 3 ;
      C = 10^12 ;
    elseif imeth == 4 ;
      C = 10^(-6) ;
    end ;  



    %  Run SVM algorithm
    %
      ker = 'linear' ;
    [nsv,alpha,bias] = SVC(mtraindat,vtrainclass,ker,C)
      %  classification version of SVM

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

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

    shpslope = -w(1) / w(2) ;
    shpyint = b / w(2) ;

    w = w / norm(w) ;




  end ;    %  of imeth if-block




  %  make next frame's graphics
  %
  correctdir = [correctdir, 0] ;
        %  correct direction vector
  mt1c = mtraindat1 * correctdir' ;
  mt2c = mtraindat2 * correctdir' ;
        %  projection of data onto correct direction

  wocorrect = w - (w * correctdir') * correctdir ;
  wocorrect = wocorrect / norm(wocorrect) ;
        %  part of w that is orthogonal to the correct direction

  mt1woc = mtraindat1 * wocorrect' ;
  mt2woc = mtraindat2 * wocorrect' ;
        %  projection of data onto discrimination direction

  wslope = w * correctdir' ;
  wslope = sqrt(1 - wslope^2) / wslope ;
  wleft = left * wslope ;
  wright = right * wslope ;


  mt1w = mtraindat1 * w' ;
  mt2w = mtraindat2 * w' ;
        %  projection of data onto discrimination direction
  h1 = bwsjpiSM(mt1w,[left; right]) ;
  h2 = bwsjpiSM(mt2w,[left; right]) ;
  h = exp(mean(log([h1; h2]))) ;
    paramstruct = struct('vxgrid',[left; right], ...
                         'vh',h) ;
  [d1kde2,kdexgrid] = kdeSM(mt1w,paramstruct) ;
  [d2kde2,kdexgrid] = kdeSM(mt2w,paramstruct) ;
  top2 = max([d1kde2; d2kde2]) ;
  bottom2 = -0.05 * top2 ;
  top2 = 1.05 * top2 ;
  d1rhts2 = bottom2 + (.7 + .1 * rhts1) * (top2 - bottom2) ;
  d2rhts2 = bottom2 + (.5 + .1 * rhts2) * (top2 - bottom2) ;
            %  random heights



  subplot(2,1,1) ;
  plot(mt1c,mt1woc,'r+',...
       mt2c,mt2woc,'bo',...
           'MarkerSize',ms,'LineWidth',lw) ;
  t1h = title([disctitstr ', dimension = ' num2str(d)],'FontSize',15) ;
  axis([left,right,bottom,top]) ;
  hold on ;
    plot([left; right],[0; 0],'m--','LineWidth',2) ;
    plot([left; right],[wleft; wright],'m-','LineWidth',2) ;
  hold off ;


  subplot(2,2,4) ;
  plot(kdexgrid,d1kde2,'r-', ...
       kdexgrid,d2kde2,'b-', ...
       [left; right],[0; 0],'m-','LineWidth',lwkde) ;
       axis([left,right,bottom2,top2]) ;
    t3h = title('Est''d Projection','FontSize',15) ;
    hold on ;
       plot(mt1w,d1rhts2,'r+','MarkerSize',ms,'LineWidth',lw) ;
       plot(mt2w,d2rhts2,'bo','MarkerSize',ms,'LineWidth',lw) ;
    hold off ;




  if iout == 1 ;    %  Increasing dimension movie

    moviestruct(d) = getframe(fighand) ;

  end ;


%pauseSM


end ;    % of iframe loop





if iout == 1 ;    %  Increasing dimension movie

  %  Reorder frames, and replay
  %
  vmorder = [(1:maxd),((maxd-1):-1:2)] ;
  moviestruct = moviestruct(vmorder) ;
          %  reorder frames, to loop back to beginning

  movie(fighand,moviestruct,nrepeat,moviefps) ;
        %  Play movie on screen


    pstr = 'DWD1figB' ;

  disp('  saving movie') ;
    movie2avi(moviestruct,pstr,'compression',moviecstr, ...
                          'keyframe',moviefps,'fps',moviefps) ;
  disp('  movie saved') ;


elseif iout == 2 ;    %  Still version for paper

  set(t1h,'FontSize',10) ;
  set(t2h,'FontSize',10) ;
  set(t3h,'FontSize',10) ;

    pstr = 'DWD1figB' ;
    orient portrait ;
    set(gcf,'PaperSize',[6.0, 4.0]) ; 
    set(gcf,'PaperPosition',[0.25, 0.25, 5.0, 3.5]) ; 
  print('-dpsc', [pstr '.ps']) ;


end ;    %  of ipart if-block





