disp('Running MATLAB script file GeomRepHDLSS1figE.m') ;
%
%    Makes Figure E for the DWD paper
%
%    Copied from HallHighDSim1.m
%
%    Which did simulations for the
%    Hall Marron High Dimensionality paper
%
%    Uses lines from 
%    and fromD:\Research\MatlabSoftware\HDD\hdd1hall1.m
%


%ipart = 2 ;     %  1 - Explore mean and variance effects
                %  2 - Compare SVM & DWD stabilities (similar to 1,
                %          but put methods on same plot)


ibw = 1 ;    %   0 for original color
             %   1 for black and white version



addpath(['C:\Documents and Settings\J S Marron\My Documents' ...
             '\Research\MatlabSoftware\HDD'],'-end') ;





  seed = 9328749087234 ;
  randn('state',seed) ;
  seed = 8459327459754 ;
  rand('state',seed) ;

  vidat = 1 ;
%  vidat = [1 2 3] ;     %  1 - N(0,I)
                        %  2 - Exp(1)
                        %  3 - 2 pt.

%  vimeth = 1 ;
%  vimeth = [1 3 7 11] ;     %  index of classification method,
  vimeth = [7 11 1 3] ;     %  index of classification method,
                           %  1 - simple mean difference
                           %  3 - 1-NN
                           %  7 - SVM, Mike Todd
                           %  11 - DWD, 1



  vd = [10 100 1000] ;
  mltypstr = strvcat(': ','--','- ') ;
  ltypechararray = char('-','--','-.',':') ;


  vmu = [0 0.5 1 1.5 3 4 5] ;
      %  Note: balance point for SVM, DWD, MD is:
      %        mu^2  =  tau^2 / n - sig^2 / m  =  20 / 16 - 4 / 16
      %              =  (20 - 4) / 16 = 1
      %        mu    =  1
      %  While balance point for NN is: 
      %        mu^2  =  tau^2 - sig^2  =  20 - 4  =  16,
      %        mu    =  4
  sig2 = 4 ;
%  vtau2 = [4 20] ;
  vtau2 = [20] ;
%  mcolstr = strvcat('b', 'g', 'm', 'r') ;
  mcolstr = strvcat('b', 'm', 'r', 'g') ;
      %  assign colors to tau's
  if ibw == 0 ;
%    vthick = [4 2 3 2] ;
    vthick = [4 2 1 1] ;
  elseif ibw == 1 ;
    vthick = [1 2 4 4] ;
  end ;

  m = 16 ;
  n = 16 ;

  ntest = 100 ;
  nsim = 1000 ;


  figure(1) ;
  clf ;
  if ibw == 1 ;
    subplot(1,3,1) ;
    subplot(1,3,2) ;
    subplot(1,3,3) ;
  end ;
  for idat = vidat ;    %  loop through data types

    if idat == 1 ;
      datstr = 'N(0,I)' ;
    elseif idat == 2 ;
      datstr = 'Exp(1)' ;
    elseif idat == 3 ;
      datstr = '2 Pt.' ;
    end ;


    it = 0 ;
    for tau2 = vtau2 ;    %  loop through tau2

      it = it + 1 ;


      %  Start plots
      %
      barht = 0.5 ;
      if ibw == 0 ;
        plot(vmu,barht * ones(length(vmu),1),'k--') ;
          left = vmu(1) ;
          right = vmu(length(vmu)) ;
          bottom = 0 ;
          topoa = 0.6 ;
          axis([left right bottom topoa]) ;
%          title(['Compare \mu, \sigma^2 = ' num2str(sig2) ...
%                        ', \tau^2 = ' num2str(tau2) ', ' datstr ...
%                        ', m = ' num2str(m) ...
%                        ', n = ' num2str(n) ...
%                        ', nsim = ' num2str(nsim) ...
%                        ', ntest = ' num2str(ntest)],'FontSize',15) ;
%          ylabel('Error Rate (overall)','FontSize',15) ;
%          xlabel('\mu','FontSize',15) ;
          ylabel('Error Rate (overall)') ;
          xlabel('\mu') ;
      end ;


      id = 0 ;
      for d = vd ;    %  loop through dimensions

        id = id + 1 ;


        %  set up graphic
        %
        if ibw == 1 ;
          subplot(1,3,id) ;
            plot(vmu,barht * ones(length(vmu),1),'k-') ;
              left = vmu(1) ;
              right = vmu(length(vmu)) ;
              bottom = 0 ;
              topoa = 0.53 ;
              axis([left right bottom topoa]) ;
              title(['d = ' num2str(d)]) ;
              if id == 1 ;
                ylabel('Error Rate (overall)') ;
              end ;
%              xlabel('\mu') ;
        end ;



        im = 0 ;
        for imeth = vimeth ;    %  loop through methods

          im = im + 1 ;


          if imeth == 1 ;
%            methstr = 'MD' ;
            methstr = 'CRD' ;
          elseif imeth == 3 ;
%            methstr = 'NN' ;
%            methstr = '1 - NN' ;
            methstr = 'NND' ;
          elseif imeth == 4 ;
            methstr = 'SVM' ;
          elseif imeth == 7 ;
            methstr = 'SVM' ;
          elseif imeth == 11 ;
            methstr = 'DWD' ;
          end ;


          vmuaperrp = [] ;
          vmuaperrn = [] ;
          vmuaperr = [] ;
          for mu = vmu ;    %  loop through mu

            disp(['   Working on mu = ' num2str(mu) ...
                      ', tau2 = ' num2str(tau2) ...
                      ', d = ' num2str(d) ...
                      ', ' methstr ', ' datstr]) ;


            vproperrp = [] ;
            vproperrn = [] ;
            vproperr = [] ;
            for isim = 1:nsim ;    %  simulation loop

              %  generate data
              %
              if idat == 1 ;
                trainp = sqrt(sig2) * randn(d,m) ;    %  the X sample
                trainn = mu + sqrt(tau2) * randn(d,n) ;    %  the Y sample
                testp = sqrt(sig2) * randn(d,ntest) ;
                testn = mu + sqrt(tau2) * randn(d,ntest) ;
              elseif idat == 2 ;
                trainp = sqrt(sig2) * (-log(rand(d,m)) - 1) ;    %  the X sample
                trainn = mu + sqrt(tau2) * (-log(rand(d,n)) - 1) ;    %  the Y sample
                testp = sqrt(sig2) * (-log(rand(d,ntest)) - 1) ;
                testn = mu + sqrt(tau2) * (-log(rand(d,ntest)) - 1) ;
              elseif idat == 3 ;
                trainp = sqrt(sig2) * (2 * (rand(d,m) > 0.5) - 1) ;
                trainn = mu + sqrt(tau2) * (2 * (rand(d,n) > 0.5) - 1) ;
                testp = sqrt(sig2) * (2 * (rand(d,ntest) > 0.5) - 1) ;
                testn = mu + sqrt(tau2) * (2 * (rand(d,ntest) > 0.5) - 1) ;
              end ;

              data = [testp testn] ;


              % do discrimination
              %
              randstate = rand('state') ;
              randnstate = randn('state') ;
                  %  save current states (since the following may reset it)

              dr = hdd1SM(data,trainp,trainn,imeth) ;

              if isempty(dr) ;
%                numerr = 2 * ntest ;
              else ;          
%                verr = (dr ~= [ones(1,ntest) -ones(1,ntest)]) ;
%                    %  one where got wrong answer
%                numerr = sum(verr) ;
                verrp = (dr(1:ntest) ~= 1) ;
                    %  one where got wrong answer for plus class
                verrn = (dr(ntest+1:2*ntest) ~= -1) ;
                    %  one where got wrong answer for minus class
                numerrp = sum(verrp) ;
                numerrn = sum(verrn) ;
                numerr = numerrp + numerrn ;
              end ;
              properrp = numerrp / ntest ;
              properrn = numerrn / ntest ;
              properr = numerr / (2 * ntest) ;
                  %  proportion of errors (on fresh data)
              vproperrp = [vproperrp; properrp] ;
              vproperrn = [vproperrn; properrn] ;
              vproperr = [vproperr; properr] ;



              rand('state',randstate) ;
              randn('state',randnstate) ;
                  %  set these back to earlier values, to get correct simulation data


            end ;    %  of isim loop


            avgproperrp = mean(vproperrp) ;
            avgproperrn = mean(vproperrn) ;
            avgproperr = mean(vproperr) ;
            vmuaperrp = [vmuaperrp; avgproperrp] ;
            vmuaperrn = [vmuaperrn; avgproperrn] ;
            vmuaperr = [vmuaperr; avgproperr] ;


          end ;    %  of mu loop


          %  Add curves to graph
          %
          figure(1) ;
          hold on ;
            if ibw == 0 ;
              plot(vmu,vmuaperr,[mcolstr(im) mltypstr(id,:)],'LineWidth',vthick(im)) ;
%              text(left + 0.2 * (right - left), ...
%                   bottom + (1 - 0.15 * im) * (topoa - bottom), ...
%                   methstr,'Color',mcolstr(im),'FontSize',15) ;
              text(left + 0.2 * (right - left), ...
                   bottom + (1 - 0.15 * im) * (topoa - bottom), ...
                   methstr,'Color',mcolstr(im)) ;
            elseif ibw == 1 ;
              plot(vmu,vmuaperr,['k' ltypechararray(im,:)],'LineWidth',vthick(im)) ;
              if id == 1
                text(left + 0.3 * (right - left), ...
                     bottom + (1 - 0.15 * im) * (topoa - bottom), ...
                     methstr,'Color','k') ;
                plot([left + 0.6 * (right - left); left + 0.9 * (right - left)], ...
                     [bottom + (1 - 0.15 * im) * (topoa - bottom); ...
                          bottom + (1 - 0.15 * im) * (topoa - bottom)], ...
                     ['k' ltypechararray(im,:)],'LineWidth',vthick(im)) ;
              end ;
            end ;
          hold off ;


        end ;    %  of imeth loop


      end ;    %  of d loop


      if ibw == 0 ;
        %  add text, etc.
        %
        figure(1) ;
        for id = 1:length(vd) ;
          hold on ;
            plot([left + 0.55 * (right - left) ; ...
                  left + 0.75 * (right - left)], ...
                 ones(2,1) * bottom + (1 - 0.15 * id) * (topoa - bottom), ...
                 ['k' mltypstr(id,:)],'LineWidth',2) ;
          hold off ;
%          text(left + 0.8 * (right - left), ...
%               bottom + (1 - 0.15 * id) * (topoa - bottom), ...
%               ['d = ' num2str(vd(id))], ...
%               'Color','k','FontSize',15) ;
          text(left + 0.8 * (right - left), ...
               bottom + (1 - 0.15 * id) * (topoa - bottom), ...
               ['d = ' num2str(vd(id))], ...
               'Color','k') ;
        end ;
      end ;


%      %  print output page
%      %
%      figure(1) ;
%      orient landscape ;
%      print('-dpsc', ['Outputs\HallHighDSim1p2t' num2str(it) ...
%                           'd' num2str(idat) '.ps']) ;


  if ibw == 0 ;
      pstr = 'GeomRepHDLSS1figE' ;
      orient portrait ;
      set(gcf,'PaperSize',[6.0, 4.1]) ; 
      set(gcf,'PaperPosition',[0.25, 0.25, 5.0, 3.6]) ; 
    print('-dpsc', [pstr '.ps']) ;
  elseif ibw == 1 ;
      pstr = 'GeomRepHDLSS1figEbw' ;
      orient portrait ;
      set(gcf,'PaperSize',[8.0, 3.1]) ; 
      set(gcf,'PaperPosition',[0.25, 0.25, 7.5, 2.6]) ; 
    print('-dps2', [pstr '.ps']) ;
  end ;




    end ;    %  of tau2 loop



  end ;    %  of idat loop



