disp('Running MATLAB script file GeomRepHDLSS1figC.m') ;
%
%    Which Makes Figure C for the 
%    "Geometric Representation of HDLSS" paper
%
%    Copied from HallHighDEg1.m
%
%    Which made Toy Examples, illustrating ideas
%    in the Hall Marron High Dimensionality paper
%


ipart = 11 ;      %  1 - n=2, N(0,I), 2-d proj.
                  %  2 - n=2, Exp(1), 2-d proj.
                  %  3 - n=2, 2 pt., 2-d proj.
                  %  11 - n=3, N(0,I), 2-d proj.
                  %  12 - n=3, Exp(1), 2-d proj.
                  %  13 - n=3, 2 pt., 2-d proj.


%ioutput = 1 ;      %  1 - Make figure for paper  (saved for later)
%                   %  2 - Make Big Pic for talk



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



vd = [2 20 200 20000] ;
%vd = [2 20 200 2000] ;
%vd = [2 10 50 250] ;
                %  Vector of Dimensions
nsim = 10 ;



if ipart == 1 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geometric Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 2, N(0,I), d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n2pd2N01' ;
  
elseif ipart == 2 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geometric Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 2, Exp(1), d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n2pd2Exp' ;
  
elseif ipart == 3 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geometric Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 2, 2 pt., d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n2pd2Pt2' ;
  
elseif ipart == 11 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geom''c Simplex Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 3, N(0,I), d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n3pd2N01' ;

elseif ipart == 12 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geom''c Simplex Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 3, Exp(1), d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n3pd2Exp' ;

elseif ipart == 13 ;

  titcellstr = {['High Dim''l Toy Data, d = ' num2str(vd(1))] ... 
                ['Geom''c Simplex Proj''n, d = ' num2str(vd(2))] ... 
                ['n = 3, 2 pt., d = ' num2str(vd(3))] ... 
                [num2str(nsim) ' replications, d = ' num2str(vd(4))]};
  pstr = 'HallHighDEg1n3pd2Pt2' ;

end ;

figure(1) ;
clf ;

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



%  Loop through dimensions
%
iplot = 0 ;
for d = vd ;

  iplot = iplot + 1 ;

  mdataplotx = [] ;
  mdataploty = [] ;
      %  matrices of data to plot

  %  simulation loop
  %
  for isim = 1:nsim ;

    if  ipart == 1  | ...
        ipart == 2  | ...
        ipart == 3  ;

      %  Generate data
      %
      if ipart == 1 ;
        mdata = randn(d,2) ;
            %  N(0,1)
      elseif ipart == 2 ;
        mdata = -log(rand(d,2)) ;
            %  Exp(1)
      elseif ipart == 3 ;
        mdata = 2 * (rand(d,2) > 0.5) - 1 ;
            %  2 point, -1,1, w.p. 0.5
      end ;

      length1 = norm(mdata(:,1)) ;
      length2 = norm(mdata(:,2)) ;
      angbtw = acos((mdata(:,1) / length1)' * (mdata(:,2) / length2)) ;
          %  lengths and angle between vectors in d space

      mproj = [length1 * [cos(pi/4 - angbtw/2); sin(pi/4 - angbtw/2)], ...
               length2 * [cos(pi/4 + angbtw/2); sin(pi/4 + angbtw/2)]] ;
          %  matrix of data projected into 2-d Geometric space
          %    1st rotated halfway back from pi/4
          %    2nd rotated halfway forward from pi/4


    elseif  ipart == 11  | ...
            ipart == 12  | ...
            ipart == 13  ;

      %  Generate data
      %
      if ipart == 11 ;
        mdata = randn(d,3) ;
            %  N(0,1)
      elseif ipart == 12 ;
        mdata = -log(rand(d,3)) ;
            %  Exp(1)
      elseif ipart == 13 ;
        mdata = 2 * (rand(d,3) > 0.5) - 1 ;
            %  2 point, -1,1, w.p. 0.5
      end ;

      vmean = mean(mdata,2) ;
      mresid = mdata - vec2matSM(vmean,3) ;

      paramstruct = struct('npc',2,...
                           'iscreenwrite',0,...
                           'viout',[0 0 0 0 1]) ;
      outstruct = pcaSM(mresid,paramstruct) ;
      mpc = getfield(outstruct,'mpc') ;
      if size(mpc,1) == 1 ;
        mpc = [mpc; [0 0 0]] ;
      end ;
           %  2x3 matrix "principal components" (scores)
           %  i.e. projections onto 2-d space

      length1 = norm(mpc(:,1)) ;
      length2 = norm(mpc(:,2)) ;
      angbtw = acos((mpc(:,1) / length1)' * (mpc(:,2) / length2)) ;
          %  lengths and angle between vectors in d space

      vmean12 = mean(mpc(:, 1:2),2) ;
          %  mean vector of first 2
      ang12 = atan2(vmean12(2),vmean12(1)) ;
          %  angle of mean12 in this plane
          %    recall atan2 want y as first argument, x as second

      angrot = 3*pi/2 - ang12 - (pi/3 - angbtw/2) ;
          %  angle to rotate, to put mean between 1 and 2 at bottom
          %  with variation in angle added, to move top point off of vertical
      mrot = [[cos(angrot) -sin(angrot)]; [sin(angrot) cos(angrot)]] ;
      mproj = mrot * mpc ;
          %  rotated version of projected data

    end ;    % ipart if-block


    mdataplotx = [mdataplotx, mproj(1,:)'] ;
    mdataploty = [mdataploty, mproj(2,:)'] ;


  end ;    % isim loop



  vax2 = axisSM([mdataplotx, mdataploty]) ;
      %  concatenate to get single value for symmetry of axes
  mindat = vax2(1) ;
  maxdat = vax2(2) ;


  %  Make plot
  %
  subplot(2,2,iplot) ;
    if ibw == 0 ;
      plot(mdataplotx,mdataploty,'o') ;
    elseif ibw == 1 ;
      plot(mdataplotx(:,1),mdataploty(:,1),'ko') ;
      hold on ;
        plot(mdataplotx(:,2),mdataploty(:,2),'kx') ;
        plot(mdataplotx(:,3),mdataploty(:,3),'k+') ;
        plot(mdataplotx(:,4),mdataploty(:,4),'k*') ;
        plot(mdataplotx(:,5),mdataploty(:,5),'ks') ;
        plot(mdataplotx(:,6),mdataploty(:,6),'kd') ;
        plot(mdataplotx(:,7),mdataploty(:,7),'kv') ;
        plot(mdataplotx(:,8),mdataploty(:,8),'k^') ;
        plot(mdataplotx(:,9),mdataploty(:,9),'k>') ;
        plot(mdataplotx(:,10),mdataploty(:,10),'k<') ;
      hold off ;
    end ;
      axis([mindat,maxdat,mindat,maxdat]) ;
      axis square ;
      hold on ;
        plot([mindat; maxdat],[0; 0],'k-') ;
        plot([0; 0],[mindat; maxdat],'k-') ;
        if  ipart ==  1  | ...
            ipart ==  2  | ...
            ipart ==  3  ;
%          plot([mindat; maxdat],[mindat; maxdat],'k--') ;
          plot([0; sqrt(d)],[sqrt(d); 0],'k--') ;
        elseif  ipart == 11  | ...
                ipart == 12  | ...
                ipart == 13  ;
          plot([0; sqrt(2*d)/2],[sqrt(6*d)/3; -sqrt(6*d)/6],'k--') ;
          plot([sqrt(2*d)/2; -sqrt(2*d)/2],[-sqrt(6*d)/6; -sqrt(6*d)/6],'k--') ;
          plot([-sqrt(2*d)/2; 0],[-sqrt(6*d)/6; sqrt(6*d)/3],'k--') ;
        end ;
      hold off ;
%      title(titcellstr{iplot},'FontSize',15) ;
      title(['d = ' num2str(d)]) ;


end ;    %  of loop through dimensions




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

elseif ibw == 1 ;
  pstr = 'GeomRepHDLSS1figCbw' ;
  orient portrait ;
  set(gcf,'PaperSize',[6.0, 5.0]) ; 
  set(gcf,'PaperPosition',[0.25, 0.25, 5.5, 4.5]) ; 
  print('-dps2', [pstr '.ps']) ;

end ;






