disp('Running MATLAB script file DWD1figD.m') ;
%
%    Makes Figure D 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 = 1 ;      %  1 - Movie of increasing dimension
                %  2 - Figure for paper



%imeth = 2 ;
     %  No imeth here (as in DWD1FigB.m & DWD1FigC.m),
     %  since DWD was not in SVMeg3.m
     %  instead use hdd1SM.m



addpath '\Research\MatlabSoftware\HDD' -end ;   
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 ;



disctitstr = 'Distance Weighted Disc.' ;
outmethstr = 'DWD' ;



%  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
  %
  imethin = 11 ;
      %  Distance Weighted Discrimination
  [dr, w] = hdd1SM(zeros(d,1),mtraindat1',mtraindat2',imethin) ;
      %  hdd1SM wants columns as data vectors
  w = w' ;
      %  return to row vector form




  %  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 = 'DWD1figD' ;

  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 = 'DWD1figD' ;
    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





