disp('Running MATLAB script file RCQ1figE.m') ;
%
%    Makes Figure E for the RCQ1 paper
%
%
%    For simluations about
%    High Dimension Discrimination
%    This explores:  uneven sample sizes
%                    different variances
%                    different dist. shapes
%                    outliers

ipart = 3 ;    %  part of program to run
               %  1 - run main simulations, as in: hdd1sim1.m
               %  2 - put things together into simulations, as in: hdd1sum1.m
               %  3 - make projection plots for 4 realizations

vidat = 14 ;
%vidat = 1:16 ;
%vidat = [3 4 7 8 11 12 15 16] ;
                  %  Vector of indices of data type
                  %  1 - same n, same var, same shape, no outliers
                  %  2 - same n, same var, same shape, outliers
                  %  3 - same n, same var, dif. shape, no outliers
                  %  4 - same n, same var, dif. shape, outliers
                  %  5 - same n, dif. var, same shape, no outliers
                  %  6 - same n, dif. var, same shape, outliers
                  %  7 - same n, dif. var, dif. shape, no outliers
                  %  8 - same n, dif. var, dif. shape, outliers
                  %  9 - same n, same var, same shape, no outliers
                  %  10 - same n, same var, same shape, outliers
                  %  11 - same n, same var, dif. shape, no outliers
                  %  12 - same n, same var, dif. shape, outliers
                  %  13 - same n, dif. var, same shape, no outliers
                  %  14 - same n, dif. var, same shape, outliers
                  %  15 - same n, dif. var, dif. shape, no outliers
                  %  16 - same n, dif. var, dif. shape, outliers
                  


% For ipart == 1  and for  ipart == 3 
%
vimeth = [32] ;
%vimeth = [1 32 111] ;
%vimeth = [32, 31] ;
%vimeth = [11 5 6] ;
%vimeth = [1 11 31 32 111]
                    %  Vector of indices on discrimination methods
                    %  1 - Mean Difference Discrimination
                    %  2 - Correlation Discrimination
                    %  5 - FLD, classical
                    %  6 - FLD, global cov.
                    %  11 - DWD 1, penalty factor = 100
                    %  12 - DWD 1, penalty factor = 10^6
                    %  21 - OSP 
                    %  31 - Robust Mean Difference
                    %  32 - Robust Mean Difference, Quantile Cut
                    %  41 - Aggregated Sign Discrimination
                    %  42 - Random Rotation, Aggregated Sign Discrimination
                    %  101 - NN - 1
                    %  102 - NN - 5
                    %  103 - NN - 11
                    %  104 - NN - 21
                    %  111 - SVM, C = 10^3
                    %  112 - SVM, C = 10^(-6)
                    %  113 - SVM, C = 10^(12)
                    %  114 - SVM, penalty factor = 100

%  For ipart == 2 
%
isum = 3 ;
                %  Vector of indices for summarization type
                %  1 - Compare MD, RMQ, SVM
                %  2 - Compare MD, DWD, RMQ, SVM
                %  3 - Compare MD, DWD, RMQ, SVM, RMD



%  For ipart == 3
%
d = 100 ;
                %  Just focus on "intermediate" case of d = 100



nsim = 100 ;
vid = 1:5 ;
vd = [10 40 100 400 1600] ;



if ipart == 1 ;    %  run main simulations


  for imeth = vimeth ;


    if imeth == 1 ;
      methstr = 'Mean Difference' ;
      methoutstr = 'MD' ;
      imethin = 1 ;
      methpar = [] ;
    elseif imeth == 2 ;
      methstr = 'Correlation' ;
      methoutstr = 'C' ;
      imethin = 2 ;
      methpar = [] ;
    elseif imeth == 5 ;
      methstr = 'FLD, Class' ;
      methoutstr = 'FLD' ;
      imethin = 5 ;
      methpar = [] ;
    elseif imeth == 6 ;
      methstr = 'FLD, Glob C' ;
      methoutstr = 'FLG' ;
      imethin = 6 ;
      methpar = [] ;
    elseif imeth == 11 ;
      methstr = 'DWD 1' ;
      methoutstr = 'DWD1' ;
      imethin = 11 ;
      methpar = [] ;
    elseif imeth == 12 ;
      methstr = 'DWD, 10^6' ;
      methoutstr = 'DWDM' ;
      imethin = 11 ;
      methpar = 10^6 ;
    elseif imeth == 21 ;
      methstr = 'OSP 1' ;
      methoutstr = 'OSP1' ;
      methpar = [] ;
      imethin = 21 ;
    elseif imeth == 31 ;
      methstr = 'RMD' ;
      methoutstr = 'RMD' ;
      methpar = [] ;
      imethin = 31 ;
    elseif imeth == 32 ;
      methstr = 'RMQ' ;
      methoutstr = 'RMQ' ;
      methpar = [] ;
      imethin = 32 ;
    elseif imeth == 41 ;
      methstr = 'ASD' ;
      methoutstr = 'ASD' ;
      methpar = [] ;
      imethin = 41 ;
    elseif imeth == 42 ;
      methstr = 'RAS' ;
      methoutstr = 'RAS' ;
      methpar = [] ;
      imethin = 42 ;
    elseif imeth == 101 ;
      methstr = '1 - NN' ;
      methoutstr = 'NN1' ;
      imethin = 3 ;
      methpar = 1 ;
    elseif imeth == 102 ;
      methstr = '5 - NN' ;
      methoutstr = 'NN5' ;
      imethin = 3 ;
      methpar = 5 ;
    elseif imeth == 103 ;
      methstr = '11 - NN' ;
      methoutstr = 'NN11' ;
      imethin = 3 ;
      methpar = 11 ;
    elseif imeth == 104 ;
      methstr = '21 - NN' ;
      methoutstr = 'NN21' ;
      imethin = 3 ;
      methpar = 21 ;
    elseif imeth == 111 ;
      C = 10^3 ;
      methstr = 'SVM, C = 10^3' ;
      methoutstr = 'SVM3' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 112 ;
      C = 10^(-6) ;
      methstr = 'SVM, C = 10^{-6}' ;
      methoutstr = 'SVMm6' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 113 ;
      C = 10^(12) ;
      methstr = 'SVM, C = 10^{12}' ;
      methoutstr = 'SVM12' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 114 ;
      methstr = 'SVM, PWD Adj.' ;
      methoutstr = 'SVMPA' ;
      imethin = 4 ;
      methpar = -100 ;
    end ;


    for idat = vidat ;    %    loop through data examples

      %  set up parameters
      %
      seed = 7543870734 ;
      rand('state',seed) ;
      randn('state',seed) ;
      if idat == 1 ;
        datstr = ['same n, same var, same shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 2 ;
        datstr = ['same n, same var, same shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 3 ;
        datstr = ['same n, same var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 4 ;
        datstr = ['same n, same var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 5 ;
        datstr = ['same n, dif. var, same shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 6 ;
        datstr = ['same n, dif. var, same shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 7 ;
        datstr = ['same n, dif. var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 8 ;
        datstr = ['same n, dif. var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 9 ;
        datstr = ['dif. n, same var, same shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 10 ;
        datstr = ['dif. n, same var, same shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 11 ;
        datstr = ['dif. n, same var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 12 ;
        datstr = ['dif. n, same var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 13 ;
        datstr = ['dif. n, dif. var, same shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 14 ;
        datstr = ['dif. n, dif. var, same shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 15 ;
        datstr = ['dif. n, dif. var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 16 ;
        datstr = ['dif. n, dif. var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      end ;
      ntestp = 4 * np ;
      ntestn = 4 * nn ;


      vmeanerr = [] ;
      vsderr = [] ;
      for id = vid ;    %  Loop through dimensions

        d = vd(id) ;
      
        vproperr = [] ;
        for isim = 1:nsim ;

          if floor(isim / 10) == (isim / 10) ;
            disp(['    Working on ' methstr ', ' datstr ', d = ' ...
                       num2str(d) ', isim = ' ...
                       num2str(isim) ' of ' num2str(nsim)]) ;
          end ;

          %  generate data
          %
          if  idat == 1  |  ...
              idat == 2  |  ...
              idat == 5  |  ...
              idat == 6  |  ...
              idat == 9  |  ...
              idat == 10  |  ...
              idat == 13  |  ...
              idat == 14  ;    %  then do same shape

            if  idat == 1  |  ...
                idat == 5  |  ...
                idat == 9  |  ...
                idat == 13  ;    % then no outliers

              trainp = sigp * (randn(d,np) + 3 / (2 * sqrt(d))) ;
              trainn = sign * (randn(d,nn) - 3 / (2 * sqrt(d))) ;
              testp = sigp * (randn(d,ntestp) + 3 / (2 * sqrt(d))) ;
              testn = sign * (randn(d,ntestn) - 3 / (2 * sqrt(d))) ;
              data = [testp testn] ;

            elseif  idat == 2  |  ...
                    idat == 6  |  ...
                    idat == 10  |  ...
                    idat == 14  ;    %  then do outliers

              trainp = sigp * (randn(d,np) + 3 / (2 * sqrt(d))) ;
              trainn = sign * (randn(d,nn) - 3 / (2 * sqrt(d))) ;
                  %  Basic Gaussian data
              outflagp = (rand(1,np) < 0.1) ;
              outflagn = (rand(1,nn) < 0.1) ;
                  %  flags 10% outliers
              trainp(1,outflagp) = trainp(1,outflagp) + 10 * sqrt(d) ;
              trainn(1,outflagn) = trainn(1,outflagn) + 10 * sqrt(d) ;
                  %  adds outlier in first coordinate only

              testp = sigp * (randn(d,ntestp) + 3 / (2 * sqrt(d))) ;
              testn = sign * (randn(d,ntestn) - 3 / (2 * sqrt(d))) ;
                  %  Basic Gaussian data
              outflagp = (rand(1,ntestp) < 0.1) ;
              outflagn = (rand(1,ntestn) < 0.1) ;
                  %  flags 10% outliers
              testp(1,outflagp) = testp(1,outflagp) + 10 * sqrt(d) ;
              testn(1,outflagn) = testn(1,outflagn) + 10 * sqrt(d) ;
                  %  adds outlier in first coordinate only

              data = [testp testn] ;



            end ;    %  of inner idat if-block


          elseif  idat == 3  |  ...
                  idat == 4  |  ...
                  idat == 7  |  ...
                  idat == 8  |  ...
                  idat == 11  |  ...
                  idat == 12  |  ...
                  idat == 15  |  ...
                  idat == 16  ;    %  then do different shape

            if  idat == 3  |  ...
                idat == 7  |  ...
                idat == 11  |  ...
                idat == 15  ;    % then no outliers

              ncutp = round(np / 2) ;
              ncutn = round(nn / 2) ;
              trainp = (sigp / 4) * (randn(d,ncutp) + 3 / (2 * sqrt(d))) ;
              trainp = [trainp (sigp * 4 * (randn(d,np-ncutp) + 3 / (2 * sqrt(d))))] ;
              trainn = (sign / 4) * (randn(d,ncutn) - 3 / (2 * sqrt(d))) ;
              trainn = [trainn (sign * 4 * (randn(d,nn-ncutn) - 3 / (2 * sqrt(d))))] ;

              ntestcutp = round(ntestp / 2) ;
              ntestcutn = round(ntestn / 2) ;
              testp = (sigp / 4) * (randn(d,ntestcutp) + 3 / (2 * sqrt(d))) ;
              testp = [testp (sigp * 4 * (randn(d,ntestp-ntestcutp) + 3 / (2 * sqrt(d))))] ;
              testn = (sign / 4) * (randn(d,ntestcutn) - 3 / (2 * sqrt(d))) ;
              testn = [testn (sign * 4 * (randn(d,ntestn-ntestcutn) - 3 / (2 * sqrt(d))))] ;
              data = [testp testn] ;

            elseif  idat == 4  |  ...
                    idat == 8  |  ...
                    idat == 12  |  ...
                    idat == 16  ;    %  then do outliers

              ncutp = round(np / 2) ;
              ncutn = round(nn / 2) ;
              trainp = (sigp / 4) * (randn(d,ncutp) + 3 / (2 * sqrt(d))) ;
              trainp = [trainp (sigp * 4 * (randn(d,np-ncutp) + 3 / (2 * sqrt(d))))] ;
              trainn = (sign / 4) * (randn(d,ncutn) - 3 / (2 * sqrt(d))) ;
              trainn = [trainn (sign * 4 * (randn(d,nn-ncutn) - 3 / (2 * sqrt(d))))] ;

              outflagp = (rand(1,np) < 0.1) ;
              outflagn = (rand(1,nn) < 0.1) ;
                  %  flags 10% outliers
              trainp(1,outflagp) = trainp(1,outflagp) + 10 * sqrt(d) ;
              trainn(1,outflagn) = trainn(1,outflagn) + 10 * sqrt(d) ;
                  %  adds outlier in first coordinate only

              ntestcutp = round(ntestp / 2) ;
              ntestcutn = round(ntestn / 2) ;
              testp = (sigp / 4) * (randn(d,ntestcutp) + 3 / (2 * sqrt(d))) ;
              testp = [testp (sigp * 4 * (randn(d,ntestp-ntestcutp) + 3 / (2 * sqrt(d))))] ;
              testn = (sign / 4) * (randn(d,ntestcutn) - 3 / (2 * sqrt(d))) ;
              testn = [testn (sign * 4 * (randn(d,ntestn-ntestcutn) - 3 / (2 * sqrt(d))))] ;

              outflagp = (rand(1,ntestp) < 0.1) ;
              outflagn = (rand(1,ntestn) < 0.1) ;
                  %  flags 10% outliers
              testp(1,outflagp) = testp(1,outflagp) + 10 * sqrt(d) ;
              testn(1,outflagn) = testn(1,outflagn) + 10 * sqrt(d) ;
                  %  adds outlier in first coordinate only

              data = [testp testn] ;


            end ;    %  of inner idat if-block


          end ;    %  of middle ipart if-block



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

          if isempty(methpar) ;
            [dr, dirvec] = hdd1SM(data,trainp,trainn,imethin) ;
          else ;
            [dr, dirvec] = hdd1SM(data,trainp,trainn,imethin,methpar) ;
          end ;

          randn('state',randnstate) ;
              %  set this back to earlier value, to get correct simulation data


          %  save results for summarization
          %
          if isempty(dr) ;
            numerr = ntest ;
          else ;          
            verr = (dr ~= [ones(1,ntestp) -ones(1,ntestn)]) ;
                %  one where got wrong answer
            numerr = sum(verr) ;
          end ;

          properr = numerr / (ntestp + ntestn) ;
              %  proportion of errors (on fresh data)
          vproperr = [vproperr; properr] ;


        end ;    %  of isim loop


        %  summarize results for this dimension
        %
        vmeanerr = [vmeanerr mean(vproperr)] ;
        vsderr = [vsderr std(vproperr)] ;


      end ;    %  of id loop



      %  save results to output file
      %
      save(['SimSummaries\HDD3sim' methoutstr ...
                  'dat' num2str(idat)], ...
                       'vmeanerr', ...
                       'vsderr', ...
                       'datstr', ...
                       'np', ...
                       'nn', ...
                       'sigp', ...
                       'sign') ;


    end ;    %  of idat loop



  end ;    %  of imeth loop



elseif ipart == 2 ;   %  make summary graphics


  if isum == 1 ;     %  Compare MD, RMQ, SVM
    vimeth = [1 32 111] ;
    colorchararray = char('r','m','b') ;
  elseif isum == 2 ;     %  Compare MD, DWD, RMQ, SVM
    vimeth = [1 32 111 11] ;
    colorchararray = char('r','m','b','g') ;
  elseif isum == 3 ;     %  Compare MD, DWD, RMQ, SVM, RMD
    vimeth = [1 32 111 11 31] ;
    colorchararray = char('r','m','b','g','c') ;
  end ;
                  %  Vector of indices on discrimination methods
                  %  1 - Mean Difference Discrimination
                  %  2 - Correlation Discrimination
                  %  11 - DWD 1, penalty factor = 100
                  %  12 - DWD 1, penalty factor = 10^6
                  %  21 - OSP 
                  %  31 - Robust Mean Difference
                  %  32 - Robust Mean Difference, Quantile Cut
                  %  41 - Aggregated Sign Discrimination
                  %  101 - NN - 1
                  %  102 - NN - 5
                  %  103 - NN - 11
                  %  104 - NN - 21
                  %  111 - SVM, C = 10^3
                  %  112 - SVM, C = 10^(-6)
                  %  113 - SVM, C = 10^(12)
                  %  114 - SVM, penalty factor = 100


  for idat = vidat ;    %    loop through data examples




    %  Start Main graphics
    %
    figure(1) ;
    clf ;

    plot([0.5; 3.5],[0.5; 0.5],'k-') ;
    axis([0.5,3.5,0,0.6]) ;
    xlabel('log_{10}(dimension)','FontSize',15) ;
    ylabel('Proportion Wrong','FontSize',15) ;


    methcount = 0 ;
    for imeth = vimeth ;    %    loop through discrimination methods

      methcount = methcount + 1 ;

      if imeth == 1 ;
        methstr = 'Mean Difference' ;
        methoutstr = 'MD' ;
        imethin = 1 ;
        methpar = [] ;
      elseif imeth == 2 ;
        methstr = 'Correlation' ;
        methoutstr = 'C' ;
        imethin = 2 ;
        methpar = [] ;
      elseif imeth == 5 ;
        methstr = 'FLD, Class' ;
        methoutstr = 'FLD' ;
        imethin = 5 ;
        methpar = [] ;
      elseif imeth == 6 ;
        methstr = 'FLD, Glob C' ;
        methoutstr = 'FLG' ;
        imethin = 6 ;
        methpar = [] ;
      elseif imeth == 11 ;
        methstr = 'DWD 1' ;
        methoutstr = 'DWD1' ;
        imethin = 11 ;
        methpar = [] ;
      elseif imeth == 12 ;
        methstr = 'DWD, 10^6' ;
        methoutstr = 'DWDM' ;
        imethin = 11 ;
        methpar = 10^6 ;
      elseif imeth == 21 ;
        methstr = 'OSP 1' ;
        methoutstr = 'OSP1' ;
        methpar = [] ;
        imethin = 21 ;
      elseif imeth == 31 ;
        methstr = 'RMD' ;
        methoutstr = 'RMD' ;
        methpar = [] ;
        imethin = 31 ;
      elseif imeth == 32 ;
        methstr = 'RMQ' ;
        methoutstr = 'RMQ' ;
        methpar = [] ;
        imethin = 32 ;
      elseif imeth == 41 ;
        methstr = 'ASD' ;
        methoutstr = 'ASD' ;
        methpar = [] ;
        imethin = 41 ;
      elseif imeth == 101 ;
        methstr = '1 - NN' ;
        methoutstr = 'NN1' ;
        imethin = 3 ;
        methpar = 1 ;
      elseif imeth == 102 ;
        methstr = '5 - NN' ;
        methoutstr = 'NN5' ;
        imethin = 3 ;
        methpar = 5 ;
      elseif imeth == 103 ;
        methstr = '11 - NN' ;
        methoutstr = 'NN11' ;
        imethin = 3 ;
        methpar = 11 ;
      elseif imeth == 104 ;
        methstr = '21 - NN' ;
        methoutstr = 'NN21' ;
        imethin = 3 ;
        methpar = 21 ;
      elseif imeth == 111 ;
        C = 10^3 ;
        methstr = 'SVM, C = 10^3' ;
        methoutstr = 'SVM3' ;
        imethin = 4 ;
        methpar = C ;
      elseif imeth == 112 ;
        C = 10^(-6) ;
        methstr = 'SVM, C = 10^{-6}' ;
        methoutstr = 'SVMm6' ;
        imethin = 4 ;
        methpar = C ;
      elseif imeth == 113 ;
        C = 10^(12) ;
        methstr = 'SVM, C = 10^{12}' ;
        methoutstr = 'SVM12' ;
        imethin = 4 ;
        methpar = C ;
      elseif imeth == 114 ;
        methstr = 'SVM, PWD Adj.' ;
        methoutstr = 'SVMPA' ;
        imethin = 4 ;
        methpar = -100 ;
      end ;



      load(['SimSummaries\HDD3sim' methoutstr ...
                'dat' num2str(idat)]) ;

      vcier = 1.96 * vsderr / sqrt(nsim - 1) ;



      %  Add these curves to plot
      %
      hold on ;
        plot(log10(vd),vmeanerr,[colorchararray(methcount) '-']) ;
        errorbar(log10(vd),vmeanerr,vcier,vcier, ...
                          [colorchararray(methcount) '-']) ;
          vax = axis ;
          tx = vax(1) + 0.1 * (vax(2) - vax(1)) ;
          ty = vax(3) + (1 - 0.1 * methcount) * (vax(4) - vax(3)) ;
        text(tx,ty,methstr,'Color',colorchararray(methcount), ...
                                        'FontSize',15) ;
      hold off ;


    end ;    %  of imeth loop

    title(datstr,'FontSize',15) ;



    %  Save this summarization plot
    %
    orient landscape ;
    print('-dpsc2',['SimPlots\HDD3sum1P' num2str(isum) ...
                    'dat' num2str(idat)]) ;



  end ;    %  of idat loop



elseif ipart == 3 ;   %  make projection plots

addpath \Research\MatlabSoftware\HDD -end ;

  for imeth = vimeth ;


    if imeth == 1 ;
      methstr = 'Mean Difference' ;
      methoutstr = 'MD' ;
      imethin = 1 ;
      methpar = [] ;
    elseif imeth == 2 ;
      methstr = 'Correlation' ;
      methoutstr = 'C' ;
      imethin = 2 ;
      methpar = [] ;
    elseif imeth == 5 ;
      methstr = 'FLD, Class' ;
      methoutstr = 'FLD' ;
      imethin = 5 ;
      methpar = [] ;
    elseif imeth == 6 ;
      methstr = 'FLD, Glob C' ;
      methoutstr = 'FLG' ;
      imethin = 6 ;
      methpar = [] ;
    elseif imeth == 11 ;
      methstr = 'DWD 1' ;
      methoutstr = 'DWD1' ;
      imethin = 11 ;
      methpar = [] ;
    elseif imeth == 12 ;
      methstr = 'DWD, 10^6' ;
      methoutstr = 'DWDM' ;
      imethin = 11 ;
      methpar = 10^6 ;
    elseif imeth == 21 ;
      methstr = 'OSP 1' ;
      methoutstr = 'OSP1' ;
      methpar = [] ;
      imethin = 21 ;
    elseif imeth == 31 ;
      methstr = 'RMD' ;
      methoutstr = 'RMD' ;
      methpar = [] ;
      imethin = 31 ;
    elseif imeth == 32 ;
      methstr = 'RMQ' ;
      methoutstr = 'RMQ' ;
      methpar = [] ;
      imethin = 32 ;
    elseif imeth == 41 ;
      methstr = 'ASD' ;
      methoutstr = 'ASD' ;
      methpar = [] ;
      imethin = 41 ;
    elseif imeth == 42 ;
      methstr = 'RAS' ;
      methoutstr = 'RAS' ;
      methpar = [] ;
      imethin = 42 ;
    elseif imeth == 101 ;
      methstr = '1 - NN' ;
      methoutstr = 'NN1' ;
      imethin = 3 ;
      methpar = 1 ;
    elseif imeth == 102 ;
      methstr = '5 - NN' ;
      methoutstr = 'NN5' ;
      imethin = 3 ;
      methpar = 5 ;
    elseif imeth == 103 ;
      methstr = '11 - NN' ;
      methoutstr = 'NN11' ;
      imethin = 3 ;
      methpar = 11 ;
    elseif imeth == 104 ;
      methstr = '21 - NN' ;
      methoutstr = 'NN21' ;
      imethin = 3 ;
      methpar = 21 ;
    elseif imeth == 111 ;
      C = 10^3 ;
      methstr = 'SVM, C = 10^3' ;
      methoutstr = 'SVM3' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 112 ;
      C = 10^(-6) ;
      methstr = 'SVM, C = 10^{-6}' ;
      methoutstr = 'SVMm6' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 113 ;
      C = 10^(12) ;
      methstr = 'SVM, C = 10^{12}' ;
      methoutstr = 'SVM12' ;
      imethin = 4 ;
      methpar = C ;
    elseif imeth == 114 ;
      methstr = 'SVM, PWD Adj.' ;
      methoutstr = 'SVMPA' ;
      imethin = 4 ;
      methpar = -100 ;
    end ;


    for idat = vidat ;    %    loop through data examples

      %  set up parameters
      %
      seed = 7543870734 ;
      rand('state',seed) ;
      randn('state',seed) ;
      if idat == 1 ;
        datstr = ['same n, same var, same shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 2 ;
        datstr = ['same n, same var, same shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 3 ;
        datstr = ['same n, same var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 4 ;
        datstr = ['same n, same var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 5 ;
        datstr = ['same n, dif. var, same shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 6 ;
        datstr = ['same n, dif. var, same shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 7 ;
        datstr = ['same n, dif. var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 8 ;
        datstr = ['same n, dif. var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 25 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 9 ;
        datstr = ['dif. n, same var, same shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 10 ;
        datstr = ['dif. n, same var, same shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 11 ;
        datstr = ['dif. n, same var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 12 ;
        datstr = ['dif. n, same var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 1 ;
      elseif idat == 13 ;
        datstr = ['dif. n, dif. var, same shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 14 ;
        datstr = ['dif. n, dif. var, same shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 15 ;
        datstr = ['dif. n, dif. var, dif. shape, no outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      elseif idat == 16 ;
        datstr = ['dif. n, dif. var, dif. shape, outliers'] ;
        np = 25 ;
        nn = 50 ;
        sigp = 1 ;
        sign = 4 ;
      end ;
      ntestp = 4 * np ;
      ntestn = 4 * nn ;


      vmeanerr = [] ;
      vsderr = [] ;

      d = 100 ;

      figure(1) ;
      clf ;

      for isim = 1:4 ;     %  just do first 4

        disp(['    Working on ' methstr ', ' datstr ', d = ' ...
                   num2str(d) ', isim = ' ...
                   num2str(isim) ' of ' num2str(4)]) ;

        %  generate data
        %
        if  idat == 1  |  ...
            idat == 2  |  ...
            idat == 5  |  ...
            idat == 6  |  ...
            idat == 9  |  ...
            idat == 10  |  ...
            idat == 13  |  ...
            idat == 14  ;    %  then do same shape

          if  idat == 1  |  ...
              idat == 5  |  ...
              idat == 9  |  ...
              idat == 13  ;    % then no outliers

            trainp = sigp * (randn(d,np) + 3 / (2 * sqrt(d))) ;
            trainn = sign * (randn(d,nn) - 3 / (2 * sqrt(d))) ;

          elseif  idat == 2  |  ...
                  idat == 6  |  ...
                  idat == 10  |  ...
                  idat == 14  ;    %  then do outliers

            trainp = sigp * (randn(d,np) + 3 / (2 * sqrt(d))) ;
            trainn = sign * (randn(d,nn) - 3 / (2 * sqrt(d))) ;
                %  Basic Gaussian data
            outflagp = (rand(1,np) < 0.1) ;
            outflagn = (rand(1,nn) < 0.1) ;
                %  flags 10% outliers
            trainp(1,outflagp) = trainp(1,outflagp) + 10 * sqrt(d) ;
            trainn(1,outflagn) = trainn(1,outflagn) + 10 * sqrt(d) ;
                %  adds outlier in first coordinate only


          end ;    %  of inner idat if-block


        elseif  idat == 3  |  ...
                idat == 4  |  ...
                idat == 7  |  ...
                idat == 8  |  ...
                idat == 11  |  ...
                idat == 12  |  ...
                idat == 15  |  ...
                idat == 16  ;    %  then do different shape

          if  idat == 3  |  ...
              idat == 7  |  ...
              idat == 11  |  ...
              idat == 15  ;    % then no outliers

            ncutp = round(np / 2) ;
            ncutn = round(nn / 2) ;
            trainp = (sigp / 4) * (randn(d,ncutp) + 3 / (2 * sqrt(d))) ;
            trainp = [trainp (sigp * 4 * (randn(d,np-ncutp) + 3 / (2 * sqrt(d))))] ;
            trainn = (sign / 4) * (randn(d,ncutn) + 3 / (2 * sqrt(d))) ;
            trainn = [trainn (sign * 4 * (randn(d,nn-ncutn) - 3 / (2 * sqrt(d))))] ;


          elseif  idat == 4  |  ...
                  idat == 8  |  ...
                  idat == 12  |  ...
                  idat == 16  ;    %  then do outliers

            ncutp = round(np / 2) ;
            ncutn = round(nn / 2) ;
            trainp = (sigp / 4) * (randn(d,ncutp) + 3 / (2 * sqrt(d))) ;
            trainp = [trainp (sigp * 4 * (randn(d,np-ncutp) + 3 / (2 * sqrt(d))))] ;
            trainn = (sign / 4) * (randn(d,ncutn) + 3 / (2 * sqrt(d))) ;
            trainn = [trainn (sign * 4 * (randn(d,nn-ncutn) - 3 / (2 * sqrt(d))))] ;

            outflagp = (rand(1,np) < 0.1) ;
            outflagn = (rand(1,nn) < 0.1) ;
                %  flags 10% outliers
            trainp(1,outflagp) = trainp(1,outflagp) + 10 * sqrt(d) ;
            trainn(1,outflagn) = trainn(1,outflagn) + 10 * sqrt(d) ;
                %  adds outlier in first coordinate only

          end ;    %  of inner idat if-block


        end ;    %  of middle ipart if-block
        
        data = [] ;



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

        if isempty(methpar) ;
          [dr, dirvec] = hdd1SM(zeros(d,1),trainp,trainn,imethin) ;
        else ;
          [dr, dirvec] = hdd1SM(zeros(d,1),trainp,trainn,imethin,methpar) ;
        end ;

        randn('state',randnstate) ;
            %  set this back to earlier value, to get correct simulation data


        %  Project data
        %
        vprojp = trainp' * dirvec ;
        vprojn = trainn' * dirvec ;


        %  compute kernel density estimates
        %
        h = bwsjpiSM([vprojp; vprojn]) ;
        vax = axisSM([vprojp; vprojn]) ;
          paramstruct = struct('vh',h, ...
                               'vxgrid',[vax(1) vax(2)]) ;
        [vkdep1, xgrid1] = kdeSM(vprojp,paramstruct) ;
        vkden1 = kdeSM(vprojn,paramstruct) ;

        if isim == 1 ;
          titstr = '4 Realizations of Projections' ;
        elseif isim == 2 ;
          titstr = datstr ;
        elseif isim == 3 ;
          titstr = methstr ;
        elseif isim == 4 ;
          titstr = ['np = ' num2str(np) ', nn = ' num2str(nn) ', d = ' num2str(d)] ;
        end ;


        %  make graphic of projected training data
        %
        subplot(2,2,isim) ;
          plot(xgrid1,vkdep1,'r-','LineWidth',2) ;
            axisSM(xgrid1,[vkdep1;vkden1]) ;
            vax = axis ;
%            title(titstr,'FontSize',12) ;
            title(['Realization ' num2str(isim)]) ;

          hold on ;
            plot(xgrid1,vkden1,'b-','LineWidth',2) ;

            plotbottom = vax(3) ;
            plottop = vax(4) ;
            yrand = plotbottom + (0.7 + 0.1 * rand(length(vprojp),1)) ...
                                                  * (plottop - plotbottom) ;
                   %  y coords for jittering
            plot(vprojp,yrand,'r+') ;

            yrand = plotbottom + (0.5 + 0.1 * rand(length(vprojn),1)) ...
                                                  * (plottop - plotbottom) ;
                   %  y coords for jittering
            plot(vprojn,yrand,'bo') ;

              tx = vax(1) + 0.1 * (vax(2) - vax(1)) ;
              ty = vax(3) + 0.9 * (vax(4) - vax(3)) ;
              voutstr = 'Class X' ;
%            text(tx,ty,voutstr,'Color','r','FontSize',18) ;
            text(tx,ty,voutstr,'Color','r') ;
              ty = vax(3) + 0.4 * (vax(4) - vax(3)) ;
              vnotoutstr = 'Class Y' ;
%            text(tx,ty,vnotoutstr,'Color','b','FontSize',18) ;
            text(tx,ty,vnotoutstr,'Color','b') ;

          hold off ;




      end ;    %  of isim loop



      %  Print Results
      %
%      orient landscape ;
%      print('-dpsc2',['SimPlots\HDD3sim1P' num2str(ipart) ...
%                      methoutstr 'dat' num2str(idat)]) ;

  pstr = 'RCQ1figE' ;
  orient portrait ;
  set(gcf,'PaperSize',[6.0, 5.0]) ; 
  set(gcf,'PaperPosition',[0.25, 0.25, 5.5, 4.5]) ; 
print('-dpsc', [pstr '.ps']) ;



    end ;    %  of idat loop



  end ;    %  of imeth loop



end ;    % of ipart if-block


