disp('Running MATLAB script file hdd1sim1.m') ;
%
%    For simluations about
%    High Dimension Discrimination

vitest = [15] ;
%vitest = [1 2 3 4 5] ;
%vitest = [11 12 13 14 15] ;
                %  Vector of indices of test type:
                %  1  - 100 simulations, with 100 tested, each time, d = 10
                %  2  - 100 simulations, with 100 tested, each time, d = 40
                %  3  - 100 simulations, with 100 tested, each time, d = 100
                %  4  - 100 simulations, with 100 tested, each time, d = 400
                %  5  - 100 simulations, with 100 tested, each time, d = 1600
                %  11 - 1 simulation, with view of projection, d = 10
                %  12 - 1 simulation, with view of projection, d = 40
                %  13 - 1 simulation, with view of projection, d = 100
                %  14 - 1 simulation, with view of projection, d = 400
                %  15 - 1 simulation, with view of projection, d = 1600


%vidat = [10 11] ;
vidat = [1 2] ;
%vidat = [1 3 5] ;
%vidat = [1 6 7] ;
%vidat = [6 7 8 9] ;
%vidat = [1 2 3 4 5 6] ;
%vidat = [1 2 3 4 5 6 7 8 9 10 11] ;
                  %  Vector of indices of data type
                  %  1 - mu = +-2.2, dim 1, d from vitest
                  %  2 - times 0.5, for d > 1, d from vitest
                  %  3 - times 2, for d > 1, d from vitest
                  %  4 - +-Chi Square 1, for d > 1, d from vitest
                  %  5 - +-log normal, for d > 1, d from vitest
                  %  6 - mu = +-2.2, dim 1, and for 20%:
                  %            +-100, dim 1 & +- 500, dim 2
                  %  7 - mu = +-2.2, dim 1, and for 10%:
                  %            +-0.1, dim 1 +- 100, same rand dim
                  %  8 - mu = +-5, others * 10
                  %  9 - donut, d/2, Embed 2, Rot'd
                  %  10 - mu = +-2, dim1, and 75% Class+ at 100
                  %  11 - mu = +-2, d1 (mirrored), 75% C+ at 100


vimeth = [12] ;
%vimeth = [1 114] ;
%vimeth = [101 102 103 104] ;
%vimeth = [111 112 113] ;
%vimeth = [1 2 11 21 101 102 103 104 111 112 113 114] ;
%vimeth = [1 11 114] ;
                    %  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 
                    %  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




if vimeth == 11 ;    %  then create file saving summaries 
                     %  of Mike Todd's "err0"

  err0outfileid = fopen('HDD1sim1err0log.txt','wt') ;
          %  'wt' is for "delete contents of this file and open 
          %               for writing" (with 't' for "text").

    titlstr1 = 'Output from the MATLAB Script HDD1sim1.m,   ' ;
    titlstr1 = [titlstr1,date] ;
  cntbytes = fprintf(err0outfileid,'%1s\n\n',titlstr1) ;

    titlstr2 = '     Records summaries of the variable "err0",' ;
  cntbytes = fprintf(err0outfileid,'%1s\n',titlstr2) ;

    titlstr3 = '     from Mike Todd''s sep.m, using SDPT3' ;
  cntbytes = fprintf(err0outfileid,'%1s\n\n\n\n',titlstr3) ;

end ;



for itest = vitest ;


  if  itest == 1  | ...
      itest == 2  | ...
      itest == 3  | ...
      itest == 4  | ...
      itest == 5  ;    %  then do simulation, testing with new data


    nsim = 100 ;

    if itest == 1 ;
      d = 10 ;
    elseif itest == 2 ;
      d = 40 ;
    elseif itest == 3 ;
      d = 100 ;
    elseif itest == 4 ;
      d = 400 ;
    elseif itest == 5 ;
      d = 1600 ;
    end ;


    for idat = vidat ;    %    loop through data examples


      for imeth = vimeth ;    %    loop through discrimination methods

        if imeth == 1 ;
          methstr = 'Mean Difference' ;
          methoutstr = 'MD' ;
          imethin = 1 ;
          methpar = [] ;
        elseif imeth == 2 ;
          methstr = 'Correlation' ;
          methoutstr = 'C' ;
          imethin = 2 ;
          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 == 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 ;


        figure(1) ;
        clf ;


        %  set up parameters
        %
        seed = 7543870734 ;
        rand('state',seed) ;
        randn('state',seed) ;
        np = 25 ;
        nn = 25 ;
        ntest = 200 ;
        if idat == 1 ;
          datstr = ['mu +-2.2, dim 1, d = ' num2str(d)] ;
        elseif idat == 2 ;
          datstr = ['mu +-2.2 & * 0.5, d = ' num2str(d)] ;
        elseif idat == 3 ;
          datstr = ['mu +-2.2 & * 2, d = ' num2str(d)] ;
        elseif idat == 4 ;
          datstr = ['mu +-2.2 & ChiSq1, d = ' num2str(d)] ;
        elseif idat == 5 ;
          datstr = ['mu +-2.2 & log-N, d = ' num2str(d)] ;
        elseif idat == 6 ;
          datstr = ['mu +-2.2 & 20% out, d = ' num2str(d)] ;
        elseif idat == 7 ;
          datstr = ['mu +-2.2 & 10% wobble, d = ' num2str(d)] ;
        elseif idat == 8 ;
          datstr = ['mu +-5 & * 10, d = ' num2str(d)] ;
        elseif idat == 9 ;
          datstr = ['Donut, Embed2, Rot, d = ' num2str(d)] ;
        elseif idat == 10 ;
          datstr = ['mu +-2, & 75% right, d = ' num2str(d)] ;
        elseif idat == 11 ;
          datstr = ['mu +-2, 75% R (Mir), d = ' num2str(d)] ;
        end ;



        vproperr = [] ;
        vangle = [] ;
        vncross0 = [] ;
        vncross = [] ;
        mdirvec = [] ;
        if vimeth == 11 ;
          verr0 = [] ;
        end ;
        for isim = 1:nsim ;

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

          %  generate data
          %
          if  idat == 1  |  ...
              idat == 2  |  ...
              idat == 3  |  ...
              idat == 4  |  ...
              idat == 5  |  ...
              idat == 6  |  ...
              idat == 7  ;

            trainp = randn(d,np) ;
            trainp(1,:) = trainp(1,:) + 2.2 ;
            trainn = randn(d,nn) ;
            trainn(1,:) = trainn(1,:) - 2.2 ;

            testp = randn(d,ntest / 2) ;
            testp(1,:) = testp(1,:) + 2.2 ;
            testn = randn(d,ntest / 2) ;
            testn(1,:) = testn(1,:) - 2.2 ;
            data = [testp testn] ;


            if idat == 2 ;
              trainp(2:d,:) = trainp(2:d,:) * 0.5 ;
              trainn(2:d,:) = trainn(2:d,:) * 0.5 ;
              data(2:d,:) = data(2:d,:) * 0.5 ;
            elseif idat == 3 ;
              trainp(2:d,:) = trainp(2:d,:) * 2 ;
              trainn(2:d,:) = trainn(2:d,:) * 2 ;
              data(2:d,:) = data(2:d,:) * 2 ;
            elseif idat == 4 ;
              randpm1p = (2 * (rand(d-1,np) > 0.5) - 1) ;
              randpm1n = (2 * (rand(d-1,nn) > 0.5) - 1) ;
              randpm1d = (2 * (rand(d-1,ntest) > 0.5) - 1) ;
                  %  takes value +-1 at random
              trainp(2:d,:) = randpm1p .* trainp(2:d,:).^2 / sqrt(3) ;
              trainn(2:d,:) = randpm1p .* trainn(2:d,:).^2 / sqrt(3) ;
              data(2:d,:) = randpm1d .* data(2:d,:).^2 / sqrt(3) ;
                  %  square to get Chi-Squared 1
                  %  and multiply by +-1, to get mean 0,
                  %  divide by sqrt(3), to make s.d. = 1
            elseif idat == 5 ;
              randpm1p = (2 * (rand(d-1,np) > 0.5) - 1) ;
              randpm1n = (2 * (rand(d-1,nn) > 0.5) - 1) ;
              randpm1d = (2 * (rand(d-1,ntest) > 0.5) - 1) ;
                  %  takes value +-1 at random
              trainp(2:d,:) = randpm1p .* exp(trainp(2:d,:)) ./ exp(1) ;
              trainn(2:d,:) = randpm1n .* exp(trainn(2:d,:)) ./ exp(1) ;
              data(2:d,:) = randpm1d .* exp(data(2:d,:)) ./ exp(1) ;
                  %  exponentiate to get standard lognormal
                  %  and multiply by +-1, to get mean 0
                  %  divide by exp(1), to make s.d. = 1
            elseif idat == 6 ;
              rand01p = rand(1,np) < 0.2 ;
              rand01n = rand(1,nn) < 0.2 ;
              rand01d = rand(1,ntest) < 0.2 ;
                  %  takes value 1 with prob. 0.2
              trainp(1:2,:) = trainp(1:2,:) + [100; 500] * rand01p ;
              data(1:2,1:(ntest/2)) = data(1:2,1:(ntest/2)) + ... 
                                  [100; 500] * rand01d(1:(ntest/2)) ;
              trainn(1:2,:) = trainn(1:2,:) - [100; 500] * rand01n ;
              data(1:2,((ntest/2)+1):ntest) = data(1:2,((ntest/2)+1):ntest) - ...
                                  [100; 500] * rand01d(((ntest/2)+1):ntest) ;
                  %  randomly (with prob 20%) add 100 or -100
            elseif idat == 7 ;

              rand01 = rand(1,np) < 0.1 ;
                  %  takes value 1 with prob. 0.1
              nrand = sum(rand01) ;
                  %  number of randomly chosen data points
              indrand = 1:np ;
              indrand = indrand(rand01) ;
                  %  indices of randomly chosen data
              vranddim = 1 + (d - 1) * rand(1,np) ;
                  %  random uniform, on scale [1,d] ;
              vranddim = ceil(vranddim) ;
                  %  vectors of dimensions to twiddle (from 2,...,d),
              for irand = 1:nrand ;
                trainp(1,indrand(irand)) = 0.1 ;
                trainp(vranddim(irand),indrand(irand)) = 100 ;
                trainn(1,indrand(irand)) = -0.1 ;
                trainn(vranddim(irand),indrand(irand)) = -100 ;
              end ;

              rand01d = rand(1,ntest/2) < 0.1 ;
                  %  takes value 1 with prob. 0.1
              nrandd = sum(rand01d) ;
                  %  number of randomly chosen data points
              indrandd = 1:ntest ;
              indrandd = indrandd(rand01d) ;
                  %  indices of randomly chosen data
              vranddim = 1 + (d - 1) * rand(1,ntest/2) ;
                  %  random uniform, on scale [1,d] ;
              vranddim = ceil(vranddim) ;
                  %  vectors of dimensions to twiddle (from 2,...,d),
              for irand = 1:nrand ;
                data(1,indrand(irand)) = 0.1 ;
                data(vranddim(irand),indrand(irand)) = 100 ;
                data(1,indrand(irand) + (ntest / 2)) = -0.1 ;
                data(vranddim(irand),indrand(irand) + (ntest / 2)) = -100 ;
              end ;

            end ;

          elseif idat == 8 ;

            trainp = randn(d,np) ;
            trainp(1,:) = trainp(1,:) + 5 ;
            trainn = randn(d,nn) ;
            trainn(1,:) = trainn(1,:) - 5 ;

            testp = randn(d,ntest / 2) ;
            testp(1,:) = testp(1,:) + 5 ;
            testn = randn(d,ntest / 2) ;
            testn(1,:) = testn(1,:) - 5 ;
            data = [testp testn] ;

            trainp(2:d,:) = trainp(2:d,:) * 10 ;
            trainn(2:d,:) = trainn(2:d,:) * 10 ;
            data(2:d,:) = data(2:d,:) * 10 ;

          elseif idat == 9 ;

            sigp = sqrt((1 + 2.2 * sqrt(2 / d)) / ...
                            (1 - 2.2 * sqrt(2 / d))) ;
                %  standard deviation to make positives have reasonable distance

            trainp = sigp * randn(d/2,np) ;
            trainn = randn(d/2,nn) ;
                %  d/2 for original data
            trainp = [(trainp.^2); trainp] ;
            trainn = [(trainn.^2); trainn] ;
                %  polynomial embedding
                %  with squares put first

            testp = sigp * randn(d/2,ntest/2) ;
            testn = randn(d/2,ntest/2) ;
                %  d/2 for original data
            testp = [(testp.^2); testp] ;
            testn = [(testn.^2); testn] ;
                %  polynomial embedding
                %  with squares put first
            data = [testp testn] ;

            mrot = (1 / sqrt(d / 2)) * ones(d / 2,1) ;
            mrot = [mrot ((1 - 1 / sqrt(d/2)) * eye(d/2, d/2 - 1))] ;
            mrot = [mrot zeros(d/2,d/2)] ;
            mrot = [mrot; [zeros(d/2,d/2) eye(d/2)]] ;
            mrot = mrot' ;
                %  rotation matrix which puts mean of squares in 1st coord.

            trainp = mrot * trainp ;
            trainn = mrot * trainn ;
            data = mrot * data ;

          elseif idat == 10 ;

            trainp = randn(d,np) ;
            trainp(1,:) = trainp(1,:) + 2.0 ;
            trainn = randn(d,nn) ;
            trainn(1,:) = trainn(1,:) - 2.0 ;

            testp = randn(d,ntest / 2) ;
            testp(1,:) = testp(1,:) + 2.0 ;
            testn = randn(d,ntest / 2) ;
            testn(1,:) = testn(1,:) - 2.0 ;
            data = [testp testn] ;

            rand01p = rand(1,np) < 0.75 ;
%            rand01d = rand(1,ntest/2) < 0.2 ;
            rand01d = rand(1,ntest/2) < 0.75 ;
                %  takes value 1 with prob. 0.75
            trainp(1,:) = trainp(1,:) + 100 * rand01p ;
            data(1,1:(ntest/2)) = data(1,1:(ntest/2)) + ... 
                                      100 * rand01d(1:(ntest/2)) ;
                %  randomly (with prob 75%) add 100 to dim 1

          elseif idat == 11 ;

            trainp = randn(d,np) ;
            trainp(1,:) = trainp(1,:) + 2.0 ;
            trainn = randn(1,nn) ;
            trainn = [trainn; trainp(2:d,:)] ;
            trainn(1,:) = trainn(1,:) - 2.0 ;

            testp = randn(d,ntest / 2) ;
            testp(1,:) = testp(1,:) + 2.0 ;
            testn = randn(1,ntest / 2) ;
            testn = [testn; testp(2:d,:)] ;
            testn(1,:) = testn(1,:) - 2.0 ;
            data = [testp testn] ;

            rand01p = rand(1,np) < 0.75 ;
%            rand01d = rand(1,ntest/2) < 0.2 ;
            rand01d = rand(1,ntest/2) < 0.75 ;
                %  takes value 1 with prob. 0.75
            trainp(1,:) = trainp(1,:) + 100 * rand01p ;
            data(1,1:(ntest/2)) = data(1,1:(ntest/2)) + ... 
                                      100 * rand01d(1:(ntest/2)) ;
                %  randomly (with prob 75%) add 100 to dim 1


          end ;



          %  count crossings, in 1st dimension
          %
          vncross0 = [vncross0; sum([(trainp(1,:) < 0) (trainn(1,:) > 0)])] ;
          vncross = [vncross; sum([(trainp(1,:) < max(trainn(1,:))) ...
                                   (trainn(1,:) > min(trainp(1,:)))])] ;


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

          if isempty(methpar) ;
            if vimeth == 11 ;
              [dr, dirvec, err0] = hdd1SM(data,trainp,trainn,imethin) ;
              verr0 = [verr0; err0] ;
            else ;
              [dr, dirvec] = hdd1SM(data,trainp,trainn,imethin) ;
            end ;
          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,ntest/2) -ones(1,ntest/2)]) ;
                %  one where got wrong answer
            numerr = sum(verr) ;
          end ;

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

          if isempty(dirvec) ;
            vangle = [] ;
          else ;
            vangle = [vangle; abs(180 * acos(dirvec(1)) / pi)] ;
                %  abs. val. angle with (1,0,...,0) in degrees
            mdirvec = [mdirvec dirvec] ;
          end ;


        end ;    %  of isim loop


        if vimeth == 11 ;    %  then write err0 summaries to file

          meanerr0 = mean(verr0) ;
          sderr0 = std(verr0) ;

            titlstr = ['For ' methoutstr ', ' datstr] ;
          cntbytes = fprintf(err0outfileid,'%1s\n',titlstr) ;

          forstr = ['       Mean of err0 = %15.10f,' ...
                    '       s.d. of err0 = %15.10f\n\n\n\n'] ;
          cntbytes = fprintf(err0outfileid,forstr,[meanerr0' sderr0]) ;

        end ;



        %  Create graphical summary
        %
        if std(vproperr) < 10^(-10) ;    %  then will have trouble with kde,
                                         %  since all answers are the same,
                                         %  so give a slight twiddle
          vproperr(1) = vproperr(1) - 0.00005 ;
        end ;
        subplot(2,2,1) ;
            paramstruct = struct('vxgrid',[0,1], ...
                                 'ndataoverlay',1, ...
                                 'titlestr',['Error Rate, ' methstr], ...
                                 'titlefontsize',15, ...
                                 'xlabelstr','Proportion Wrong', ...
                                 'ylabelstr','density', ...
                                 'labelfontsize',15) ;
          kdeSM(vproperr,paramstruct) ;
          hold on ;
            vax = axis ;
              tx = vax(1) + 0.10 * (vax(2) - vax(1)) ;
              ty = vax(3) + 0.70 * (vax(4) - vax(3)) ;
            text(tx,ty,['Average Error Rate = ' num2str(mean(vproperr))], ...
                           'FontSize',15) ;
          hold off ;

        if ~isempty(vangle) ;    %  then make angle plot
          subplot(2,2,3) ;
            if std(vangle) < 10^(-10) ;    %  then will have trouble with kde,
                                             %  since all answers are the same,
                                             %  so give a slight twiddle
              vangle(1) = vangle(1) + 0.0001 ;
            end ;
              paramstruct = struct('vxgrid',[0,90], ...
                                   'ndataoverlay',1, ...
                                   'titlestr',['Angle, ' datstr], ...
                                   'titlefontsize',15, ...
                                   'xlabelstr','Angle from Optimal (degrees)', ...
                                   'ylabelstr','density', ...
                                   'labelfontsize',15) ;
            kdeSM(vangle,paramstruct) ;
            hold on ;
              vax = axis ;
                tx = vax(1) + 0.10 * (vax(2) - vax(1)) ;
                ty = vax(3) + 0.70 * (vax(4) - vax(3)) ;
              text(tx,ty,['Average Angle = ' num2str(mean(vangle))], ...
                           'FontSize',15) ;
            hold off ;

        else ;
          subplot(2,2,3) ;
          title(datstr,'FontSize',15) ;
        end ;

        maxvnc0 = max([vncross0; 5]) ;
        subplot(2,2,2) ;
          hist(vncross0,[0:maxvnc0]') ;
            title('# in training data, that cross 0','FontSize',15) ;
            xlabel('Number of 0 Crossings','FontSize',15) ;
            ylabel('Counts','FontSize',15) ;
            vax = axis ;
            axis([-0.5 (maxvnc0 + 0.5) vax(3) vax(4)]) ;
            vax = axis ;
              tx = vax(1) + 0.3 * (vax(2) - vax(1)) ;
              ty = vax(3) + 0.9 * (vax(4) - vax(3)) ;
            text(tx,ty,['n_p = ' num2str(np) '      n_n = ' num2str(nn)], ...
                                      'FontSize',15) ;
              ty = vax(3) + 0.8 * (vax(4) - vax(3)) ;
            text(tx,ty,['n_{test} = ' num2str(ntest) ...
                        '       n_{sim} = ' num2str(nsim)], ...
                                      'FontSize',15) ;

        maxvnc = max([vncross; 5]) ;
          subplot(2,2,4) ;
            hist(vncross,[0:maxvnc]') ;
              title('# in training data, that overlap other class', ...
                                       'FontSize',15) ;
              xlabel('Number of Crossings','FontSize',15) ;
              ylabel('Counts','FontSize',15) ;
              vax = axis ;
              axis([-0.5 (maxvnc + 0.5) vax(3) vax(4)]) ;

          orient landscape ;
        print('-dpsc2',['SimPlots\HDD1sim' methoutstr ...
                        'dat' num2str(idat) ...
                        'd' num2str(d)]) ;


        save(['SimSummaries\HDD1sim' methoutstr ...
                  'dat' num2str(idat) 'd' num2str(d)], ...
                       'vproperr', ...
                       'vangle', ...
                       'mdirvec') ;


      end ;    %  of imeth loop


    end ;    %  of idat loop



  elseif  itest == 11  | ...
          itest == 12  | ...
          itest == 13  | ...
          itest == 14  | ...
          itest == 15  ;    %  then make projection view


    if itest == 11 ;
      d = 10 ;
    elseif itest == 12 ;
      d = 40 ;
    elseif itest == 13 ;
      d = 100 ;
    elseif itest == 14 ;
      d = 400 ;
    elseif itest == 15 ;
      d = 1600 ;
    end ;



    for idat = vidat ;    %    loop through data examples


      for imeth = vimeth ;    %    loop through discrimination methods

        if imeth == 1 ;
          methstr = 'Mean Difference' ;
          methoutstr = 'MD' ;
          imethin = 1 ;
          methpar = [] ;
        elseif imeth == 2 ;
          methstr = 'Correlation' ;
          methoutstr = 'C' ;
          imethin = 2 ;
          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 == 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 ;


        %  set up parameters
        %
        seed = 7543870734 ;
        randn('state',seed) ;
        np = 25 ;
        nn = 25 ;
        ntest = 200 ;
        if idat == 1 ;
          datstr = ['mu +-2.2, dim 1, d = ' num2str(d)] ;
        elseif idat == 2 ;
          datstr = ['mu +-2.2 & * 0.5, d = ' num2str(d)] ;
        elseif idat == 3 ;
          datstr = ['mu +-2.2 & * 2, d = ' num2str(d)] ;
        elseif idat == 4 ;
          datstr = ['mu +-2.2 & ChiSq1, d = ' num2str(d)] ;
        elseif idat == 5 ;
          datstr = ['mu +-2.2 & log-N, d = ' num2str(d)] ;
        elseif idat == 6 ;
          datstr = ['mu +-2.2 & 20% out, d = ' num2str(d)] ;
        elseif idat == 7 ;
          datstr = ['mu +-2.2 & 10% wobble, d = ' num2str(d)] ;
        elseif idat == 8 ;
          datstr = ['mu +-5 & * 10, d = ' num2str(d)] ;
        elseif idat == 9 ;
          datstr = ['Donut, Embed2, Rot, d = ' num2str(d)] ;
        elseif idat == 10 ;
          datstr = ['mu +-2, & 75% right, d = ' num2str(d)] ;
        elseif idat == 11 ;
          datstr = ['mu +-2, 75% R (Mir), d = ' num2str(d)] ;
        end ;


        figure(1) ;
        clf ;


        %  generate data
        %
        if  idat == 1  |  ...
            idat == 2  |  ...
            idat == 3  |  ...
            idat == 4  |  ...
            idat == 5  |  ...
            idat == 6  |  ...
            idat == 7  ;

          trainp = randn(d,np) ;
          trainp(1,:) = trainp(1,:) + 2.2 ;
          trainn = randn(d,nn) ;
          trainn(1,:) = trainn(1,:) - 2.2 ;

          testp = randn(d,ntest / 2) ;
          testp(1,:) = testp(1,:) + 2.2 ;
          testn = randn(d,ntest / 2) ;
          testn(1,:) = testn(1,:) - 2.2 ;
          data = [testp testn] ;


          if idat == 2 ;
            trainp(2:d,:) = trainp(2:d,:) * 0.5 ;
            trainn(2:d,:) = trainn(2:d,:) * 0.5 ;
            data(2:d,:) = data(2:d,:) * 0.5 ;
          elseif idat == 3 ;
            trainp(2:d,:) = trainp(2:d,:) * 2 ;
            trainn(2:d,:) = trainn(2:d,:) * 2 ;
            data(2:d,:) = data(2:d,:) * 2 ;
          elseif idat == 4 ;
            randpm1p = (2 * (rand(d-1,np) > 0.5) - 1) ;
            randpm1n = (2 * (rand(d-1,nn) > 0.5) - 1) ;
            randpm1d = (2 * (rand(d-1,ntest) > 0.5) - 1) ;
                %  takes value +-1 at random
            trainp(2:d,:) = randpm1p .* trainp(2:d,:).^2 / sqrt(3) ;
            trainn(2:d,:) = randpm1p .* trainn(2:d,:).^2 / sqrt(3) ;
            data(2:d,:) = randpm1d .* data(2:d,:).^2 / sqrt(3) ;
                %  square to get Chi-Squared 1
                %  and multiply by +-1, to get mean 0,
                %  divide by sqrt(3), to make s.d. = 1
          elseif idat == 5 ;
            randpm1p = (2 * (rand(d-1,np) > 0.5) - 1) ;
            randpm1n = (2 * (rand(d-1,nn) > 0.5) - 1) ;
            randpm1d = (2 * (rand(d-1,ntest) > 0.5) - 1) ;
                %  takes value +-1 at random
            trainp(2:d,:) = randpm1p .* exp(trainp(2:d,:)) ./ exp(1) ;
            trainn(2:d,:) = randpm1n .* exp(trainn(2:d,:)) ./ exp(1) ;
            data(2:d,:) = randpm1d .* exp(data(2:d,:)) ./ exp(1) ;
                %  exponentiate to get standard lognormal
                %  and multiply by +-1, to get mean 0
                %  divide by exp(1), to make s.d. = 1
          elseif idat == 6 ;
            rand01p = rand(1,np) < 0.2 ;
            rand01n = rand(1,nn) < 0.2 ;
            rand01d = rand(1,ntest) < 0.2 ;
                %  takes value 1 with prob. 0.2
            trainp(1:2,:) = trainp(1:2,:) + [100; 500] * rand01p ;
            data(1:2,1:np) = data(1:2,1:np) + ... 
                                [100; 500] * rand01d(1:np) ;
            trainn(1:2,:) = trainn(1:2,:) - [100; 500] * rand01n ;
            data(1:2,(np+1):(np+nn)) = data(1:2,(np+1):(np+nn)) - ...
                                [100; 500] * rand01d((np+1):(np+nn)) ;
                %  randomly (with prob 20%) add 100 or -100
          elseif idat == 7 ;

            rand01 = rand(1,np) < 0.1 ;
                %  takes value 1 with prob. 0.1
            nrand = sum(rand01) ;
                %  number of randomly chosen data points
            indrand = 1:np ;
            indrand = indrand(rand01) ;
                %  indices of randomly chosen data
            vranddim = 1 + (d - 1) * rand(1,np) ;
                %  random uniform, on scale [1,d] ;
            vranddim = ceil(vranddim) ;
                %  vectors of dimensions to twiddle (from 2,...,d),
            for irand = 1:nrand ;
              trainp(1,indrand(irand)) = 0.1 ;
              trainp(vranddim(irand),indrand(irand)) = 100 ;
              trainn(1,indrand(irand)) = -0.1 ;
              trainn(vranddim(irand),indrand(irand)) = -100 ;
            end ;

            rand01d = rand(1,ntest/2) < 0.1 ;
                %  takes value 1 with prob. 0.1
            nrandd = sum(rand01d) ;
                %  number of randomly chosen data points
            indrandd = 1:ntest ;
            indrandd = indrandd(rand01d) ;
                %  indices of randomly chosen data
            vranddim = 1 + (d - 1) * rand(1,ntest/2) ;
                %  random uniform, on scale [1,d] ;
            vranddim = ceil(vranddim) ;
                %  vectors of dimensions to twiddle (from 2,...,d),
            for irandd = 1:nrandd ;
              data(1,indrandd(irandd)) = 0.1 ;
              data(vranddim(irandd),indrandd(irandd)) = 100 ;
              data(1,indrandd(irandd) + (ntest / 2)) = -0.1 ;
              data(vranddim(irandd),indrandd(irandd) + (ntest / 2)) = -100 ;
            end ;

          end ;

        elseif idat == 8 ;

          trainp = randn(d,np) ;
          trainp(1,:) = trainp(1,:) + 5 ;
          trainn = randn(d,nn) ;
          trainn(1,:) = trainn(1,:) - 5 ;

          testp = randn(d,ntest / 2) ;
          testp(1,:) = testp(1,:) + 5 ;
          testn = randn(d,ntest / 2) ;
          testn(1,:) = testn(1,:) - 5 ;
          data = [testp testn] ;

          trainp(2:d,:) = trainp(2:d,:) * 10 ;
          trainn(2:d,:) = trainn(2:d,:) * 10 ;
          data(2:d,:) = data(2:d,:) * 10 ;

        elseif idat == 9 ;

          sigp = sqrt((1 + 2.2 * sqrt(2 / d)) / ...
                          (1 - 2.2 * sqrt(2 / d))) ;
              %  standard deviation to make positives have reasonable distance

          trainp = sigp * randn(d/2,np) ;
          trainn = randn(d/2,nn) ;
              %  d/2 for original data
          trainp = [(trainp.^2); trainp] ;
          trainn = [(trainn.^2); trainn] ;
              %  polynomial embedding
              %  with squares put first

          testp = sigp * randn(d/2,ntest/2) ;
          testn = randn(d/2,ntest/2) ;
              %  d/2 for original data
          testp = [(testp.^2); testp] ;
          testn = [(testn.^2); testn] ;
              %  polynomial embedding
              %  with squares put first
          data = [testp testn] ;

          mrot = (1 / sqrt(d / 2)) * ones(d / 2,1) ;
          mrot = [mrot ((1 - 1 / sqrt(d/2)) * eye(d/2, d/2 - 1))] ;
          mrot = [mrot zeros(d/2,d/2)] ;
          mrot = [mrot; [zeros(d/2,d/2) eye(d/2)]] ;
          mrot = mrot' ;
              %  rotation matrix which puts mean of squares in 1st coord.

          trainp = mrot * trainp ;
          trainn = mrot * trainn ;
          data = mrot * data ;
              %  rotated version of data

          cent = ((1 + sigp^2) / 2) * d ;
          trainp(1,:) = trainp(1,:) - cent ;
          trainn(1,:) = trainn(1,:) - cent ;
          data(1,:) = data(1,:) - cent ;
              %  recenter first coordinate

        elseif idat == 10 ;

          trainp = randn(d,np) ;
          trainp(1,:) = trainp(1,:) + 2.0 ;
          trainn = randn(d,nn) ;
          trainn(1,:) = trainn(1,:) - 2.0 ;

          testp = randn(d,ntest / 2) ;
          testp(1,:) = testp(1,:) + 2.0 ;
          testn = randn(d,ntest / 2) ;
          testn(1,:) = testn(1,:) - 2.0 ;
          data = [testp testn] ;

          rand01p = rand(1,np) < 0.75 ;
%          rand01d = rand(1,ntest/2) < 0.2 ;
          rand01d = rand(1,ntest/2) < 0.75 ;
              %  takes value 1 with prob. 0.75
          trainp(1,:) = trainp(1,:) + 100 * rand01p ;
          data(1,1:(ntest/2)) = data(1,1:(ntest/2)) + ... 
                                    100 * rand01d(1:(ntest/2)) ;
              %  randomly (with prob 75%) add 100 to dim 1

        elseif idat == 11 ;

          trainp = randn(d,np) ;
          trainp(1,:) = trainp(1,:) + 2.0 ;
          trainn = randn(1,nn) ;
          trainn = [trainn; trainp(2:d,:)] ;
          trainn(1,:) = trainn(1,:) - 2.0 ;

          testp = randn(d,ntest / 2) ;
          testp(1,:) = testp(1,:) + 2.0 ;
          testn = randn(1,ntest / 2) ;
          testn = [testn; testp(2:d,:)] ;
          testn(1,:) = testn(1,:) - 2.0 ;
          data = [testp testn] ;

          rand01p = rand(1,np) < 0.75 ;
%          rand01d = rand(1,ntest/2) < 0.2 ;
          rand01d = rand(1,ntest/2) < 0.75 ;
              %  takes value 1 with prob. 0.75
          trainp(1,:) = trainp(1,:) + 100 * rand01p ;
          data(1,1:(ntest/2)) = data(1,1:(ntest/2)) + ... 
                                    100 * rand01d(1:(ntest/2)) ;
              %  randomly (with prob 75%) add 100 to dim 1

        end ;



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


        verr = (dr ~= [ones(1,ntest/2) -ones(1,ntest/2)]) ;
            %  one where got wrong answer
        numerr = sum(verr) ;
        properr = numerr / ntest ;

        ang = abs(180 * acos(dirvec(1)) / pi) ;
                %  abs. val. angle with (1,0,...,0) in degrees


        %  make main graphic
        %
        correctdir = [1; zeros(d-1,1)] ;
              %  correct direction vector
        mt1c = trainp' * correctdir ;
        mt2c = trainn' * correctdir ;
              %  projection of data onto correct direction


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

        mt1woc = trainp' * wocorrect ;
        mt2woc = trainn' * wocorrect ;
              %  projection of data onto discrimination direction

        mt1w = trainp' * dirvec ;
        mt2w = trainn' * dirvec ;
              %  projection of data onto discrimination direction

        left = min([mt1c; mt2c; mt1woc; mt2woc; mt1w; mt2w]) ;
        right = max([mt1c; mt2c; mt1woc; mt2woc; mt1w; mt2w]) ;
        range = right - left ;
        left = left - 0.05 * range ;
        right = right + 0.05 * range ;

        bottom = left ;
        top = right ;

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


        %  Compute needed density estimates, etc.
        %
          paramstruct = struct('vxgrid',[left; right]) ;
        [d1kde1,kdexgrid] = kdeSM(mt1c,paramstruct) ;
        [d2kde1,kdexgrid] = kdeSM(mt2c,paramstruct) ;
        top1 = max([d1kde1; d2kde1]) ;
        bottom1 = -0.05 * top1 ;
        top1 = 1.05 * top1 ;
        rhts1 = rand(np,1) ;
        rhts2 = rand(nn,1) ;
        d1rhts1 = bottom1 + (.7 + .1 * rhts1) * (top1 - bottom1) ;
        d2rhts1 = bottom1 + (.5 + .1 * rhts2) * (top1 - bottom1) ;
                    %  random heights

          paramstruct = struct('vxgrid',[left; right]) ;
        [d1kde2,kdexgrid] = kdeSM(mt1w,paramstruct) ;
        [d2kde2,kdexgrid] = kdeSM(mt2w,paramstruct) ;
        top2 = max([d1kde2; d2kde2]) ;
        bottom2 = -0.05 * top2 ;
        top2 = 1.05 * top2 ;
        rhts1 = rand(np,1) ;
        rhts2 = rand(nn,1) ;
        d1rhts2 = bottom2 + (.7 + .1 * rhts1) * (top2 - bottom2) ;
        d2rhts2 = bottom2 + (.5 + .1 * rhts2) * (top2 - bottom2) ;
                    %  random heights


        %  do actual plots
        %
        subplot(2,1,1) ;
        plot(mt1c,mt1woc,'r+',...
             mt2c,mt2woc,'bo',...
                 'MarkerSize',10,'LineWidth',3) ;
        title([methstr ',    ' datstr ',    2-d Projection'],'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) ;
            tx = left + 0.1 * (right - left) ;
            ty = bottom + 0.85 * (top - bottom) ;
          text(tx,ty,['Error Rate = ' num2str(properr) ...
                      '    (based on ' num2str(ntest) ...
                      ' new data),    Abs. Angle = ' ...
                      num2str(ang,3) ' (deg)'], ...
                      'FontSize',15) ;
        hold off ;


        subplot(2,2,3) ;
        plot(kdexgrid,d1kde1,'r-', ...
             kdexgrid,d2kde1,'b-', ...
             [left; right],[0; 0],'m--','LineWidth',2) ;
             axis([left,right,bottom1,top1]) ;
          title('True Projection','FontSize',15) ;
          hold on ;
             plot(mt1c,d1rhts1,'r+','MarkerSize',10,'LineWidth',3) ;
             plot(mt2c,d2rhts1,'bo','MarkerSize',10,'LineWidth',3) ;
          hold off ;


        subplot(2,2,4) ;
        plot(kdexgrid,d1kde2,'r-', ...
             kdexgrid,d2kde2,'b-', ...
             [left; right],[0; 0],'m-','LineWidth',2) ;
             axis([left,right,bottom2,top2]) ;
          title('Projection on Disc''n Direction Vector','FontSize',15) ;
          hold on ;
             plot(mt1w,d1rhts2,'r+','MarkerSize',10,'LineWidth',3) ;
             plot(mt2w,d2rhts2,'bo','MarkerSize',10,'LineWidth',3) ;
          hold off ;


          orient landscape ;
        print('-dpsc2',['SimDir\HDD1DirPlot' methoutstr ...
                        'dat' num2str(idat) ...
                        'd' num2str(d)]) ;



      end ;    %  of imeth loop


    end ;    %  of idat loop


  end ;    %  of itest if-block


end ;    %  of itest loop



if vimeth == 11 ;    %  then did file saving summaries 
                     %  of Mike Todd's "err0"

  fclose(err0outfileid) ;
      %  !!!! Careful !!!!  Need to do this to be able
      %  to work with files after running matlab. 

end ;


