function makeplot = ViewDimAway(mdata,vtarget,paramstruct) 
% VIEWDIMAWAY, View of DIMensions AWAY
%    Studies "range" of high dimensional point clouds,
%    by graphics which focus on "dimensionality" of 
%    spherical normal point clouds.
%    Looks at "angles" on the scale of "dimensions away", and
%    compares distributions to those expected assuming a
%    lower dimensional spherical normal distribution.
%    Useful for studying "stability" of high dimensional
%    statistical methods.
%    To only calculate "dimensions away", use VecDimAway.m
% Inputs:
%   mdata       - d x n matrix of data,
%                      organized as n column vectors, 
%                          of d-dim data vectors
%
%   vtarget     - target vector, work with angles between
%                    this and mdata.  Use pairwise angles when 
%                    this is set to 0 or not specified
%
%   paramstruct - a Matlab structure of input parameters
%                    Use: "help struct" and "help datatypes" to
%                         learn about these.
%                    Create one, using commands of the form:
%
%       paramstruct = struct('field1',values1,...
%                            'field2',values2,...
%                                             ) ;
%
%                          where any of the following can be used,
%                          these are optional, misspecified values
%                          revert to defaults
%
%    fields            values
%
%    iout             1  (default)  use 3 panels:  DimAway and kde,
%                                      Prob. Int. transformed,
%                                      and Gaussian Q-Q plot
%                     2  use 2 panels: DimAway and P. I. T.
%                     3  DimAway and kde only
%                     4  Prob. Int. transformed only
%                     5  Gaussian Q-Q plot only
%
%    imovie           0  (default)  make a single still plot
%                            (reset to this for iout > 2)
%                     1  make movie version
%
%    d0               dimension to work with:
%                       for imovie = 0,  use this value
%                            when this is 0, estimate from data
%                            using sample standard deviation
%                               (default is 0, i.e. estimate)
%                       for imovie = 1,  use d0,...,d
%                            (default is 1)
%
%    dstart          (for movie version only) start of d range
%    dend            corresponding end of d range
%                            (these override d0 above when specified)
%
%    savestr          string controlling saving of output,
%                         either a full path, or a file prefix to
%                         save in matlab's current directory
%                     unspecified:  results only appear on screen
%                     movie version (imovie = 1):
%                         add .mpg and create MPEG file
%                     static version (imovie = 0):  
%                         add .ps, and save as color postscripT
%
%    ndataoverlay     0  no data plot
%                     1  (default) plot 1000 points (random choice)
%                     2  plot full data set
%                     n > 2  plot n random points
%
%    iscreenwrite     0  (default)  no screen writes
%                     1  write to screen to show progress
%
%    nrepeat          number of times to repeat movie (default = 5)
%
%    titlestr         Title String
%                         iout = 1,2: add to "Dim. Away" in top plot
%                         iout = 3: add to "Dim. Away"
%                         iout = 4: add to "d0-d Gauss PIT"
%                         iout = 5: add to "d0-d Gauss Q-Q"
%
%    vfigpos          Screen Coordinates for Figure Position
%                     Could try: [600 150 400 550]  for iout = 1
%                                [550 300 450 400]  for iout = 2
%                     Default: use current Figure Position setting
%
%
% Outputs:
%     For iout = 1,2:   graphics in current Figure
%     For iout = 3,4,5:   graphics in current axes
%     When savestr exists,
%        For imovie = 1:  MPEG file saved in 'savestr'.mpg
%        For imovie = 0:  Postscript file saved in 'savestr'.ps
%
% Assumes path can find personal functions:
%    VecDimAway.m
%    bwsjpib.m
%    gpkde.m

%    Copyright (c) J. S. Marron 1999



%  First check whether have target vector, 
%  or are doing pairwise analysis
if nargin >= 2 ;
  if vtarget == 0 ;
    itarget = 0 ;
          %  do pairwise analysis
  else ;
    itarget = 1 ;
          %  do analysis comparing data to vtarget
  end ;
else ;
  itarget = 0 ;
          %  do pairwise analysis
end ;




d = size(mdata,1) ;
          %  dimensionality of data vectors (# of rows)
n = size(mdata,2) ;
          %  number of data vectors (# of cols)




%  Set all parameters to defaults
%
iout = 1 ;
imovie = 0 ;
d0 = 0 ;
%  note no defaults for dstart and dend
savestr = [] ;
ndataoverlay = 1 ;
iscreenwrite = 0 ;
nrepeat = 5 ;
titlestr = [] ;
vfigpos = [] ;





%  Now update parameters as specified,
%  by parameter structure (if it is used)
%
if nargin >=3 ;   %  then paramstruct has been added

  if isfield(paramstruct,'iout') ;    %  then change to input value
    iout = getfield(paramstruct,'iout') ; 
  end ;

  if isfield(paramstruct,'imovie') ;    %  then change to input value
    imovie = getfield(paramstruct,'imovie') ; 
  end ;

  if imovie == 1 ;  %  then change default d0
    d0 = 1 ;
  end ;

  if isfield(paramstruct,'d0') ;    %  then change to input value
    d0 = getfield(paramstruct,'d0') ; 
  end ;

  if isfield(paramstruct,'dstart') ;    %  then change to input value
    dstart = getfield(paramstruct,'dstart') ; 
    if imovie == 1 ;
      d0 = dstart ;
    end ;
  else ;    %  need to set to default value
    dstart = d0 ;
  end ;

  if isfield(paramstruct,'dend') ;    %  then change to input value
    dend = getfield(paramstruct,'dend') ; 
  else ;    %  need to set to default value
    dend = d ;
  end ;

  if isfield(paramstruct,'savestr') ;    %  then use input value
    savestr = getfield(paramstruct,'savestr') ; 
    if ~ischar(savestr) ;    %  then invalid input, so give warning
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from ViewDimAway.m:     !!!') ;
      disp('!!!   Invalid savestr,           !!!') ;
      disp('!!!   using default of no save   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      savestr = [] ;
    end ;
  end ;

  if isfield(paramstruct,'ndataoverlay') ;    %  then change to input value
    ndataoverlay = getfield(paramstruct,'ndataoverlay') ; 
  end ;

  if isfield(paramstruct,'iscreenwrite') ;    %  then change to input value
    iscreenwrite = getfield(paramstruct,'iscreenwrite') ; 
  end ;

  if isfield(paramstruct,'nrepeat') ;    %  then change to input value
    nrepeat = getfield(paramstruct,'nrepeat') ; 
  end ;

  if isfield(paramstruct,'titlestr') ;    %  then change to input value
    titlestr = getfield(paramstruct,'titlestr') ; 
  end ;

  if isfield(paramstruct,'vfigpos') ;    %  then change to input value
    vfigpos = getfield(paramstruct,'vfigpos') ; 
  end ;


end ;  %  of resetting of input parameters




%  Set other parameters
%
seed = 43845633 ;
rand('seed',seed) ;
          %  used by jitter plots




%  Turn iout into control parameters
%
if iout == 1 ;
  viplot = [1; 1; 1] ;
         %  indicators for ordinary, transformed, Q-Q
elseif iout == 2 ;
  viplot = [1; 1; 0] ;
elseif iout == 3 ;
  viplot = [1; 0; 0] ;
elseif iout == 4 ;
  viplot = [0; 1; 0] ;
elseif iout == 5 ;
  viplot = [0; 0; 1] ;
else ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Warning from ViewDimAway.m:      !!!') ;
  disp('!!!   bad value of iout                !!!') ;
  disp('!!!   resetting to iout = 1 (default)  !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  viplot = [1; 1; 1] ;
end ;
nplot = sum(viplot) ;


if  nplot == 1  &  imovie == 1  ;
          %  then reset imovie, since it focuses on figures, not axes 
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Warning from ViewDimAway.m:   !!!') ;
  disp('!!!   resetting imovie to 0         !!!') ;
  disp('!!!   (use iout <= 2 for movie)     !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  imovie = 0 ;
end ;





%  set movie parameters
%
nrif = 4 ;
          %  number of times to Repeat Individual Movie Frames
          %    (controls speed of .mpeg version)
movopt = 1 ;
%  Movie options, currently using same as Cornean
%
%  This is MATLAB "help" to use for choosing parameters:
%	MPGWRITE(M, map, 'filename', options) Encodes M in MPEG format
% 	using the specified colormap and writes the result to the
% 	specified file.  The options argument is an optional vector of
% 	8 or fewer options where each value has the following meaning:
% 	1. REPEAT:
% 		An integer number of times to repeat the movie
% 		(default is 1).
% 	2. P-SEARCH ALGORITHM:
% 		0 = logarithmic	(fastest, default value)
% 		1 = subsample
% 		2 = exhaustive	(better, but slow)
% 	3. B-SEARCH ALGORITHM:
% 		0 = simple	(fastest)
% 		1 = cross2	(slightly slower, default value)
% 		2 = exhaustive	(very slow)
% 	4. REFERENCE FRAME:
% 		0 = original	(faster, default)
% 		1 = decoded	(slower, but results in better quality)
% 	5. RANGE IN PIXELS:
% 		An integer search radius.  Default is 10.
% 	6. I-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 8.
% 	7. P-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 10.
% 	8. B-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 25.





%  check input data OK
%
if n < 2 ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from ViewDimAway.m:      !!!') ;
  disp('!!!   Need more than 1 col in mdata  !!!') ;
  disp('!!!   Terminating Execution          !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  return ;
end ;

if itarget == 1 ;    %  then have a target vector
  if size(vtarget,2) > 1 ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from ViewDimAway.m:      !!!') ;
    disp('!!!   Need vtarget a col vector      !!!') ;
    disp('!!!   Terminating Execution          !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    return ;
  end ;

  if length(vtarget) ~= d ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from ViewDimAway.m:       !!!') ;
    disp('!!!   vtarget needs same dim as data  !!!') ;
    disp('!!!   Terminating Execution           !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    return ;
  end ;
end ;





%  Do Dimensions Away Calculation
%
if itarget == 1 ;
  dimaway = VecDimAway(mdata,vtarget,1) ;
          %  1 to do normalization to unit vectors operation
  dimaway = dimaway' ;
else ;
  dimaway = VecDimAway(mdata,1) ;
          %  1 to do normalization to unit vectors operation

  flag = logical(1-triu(ones(n))) ;
          %  matrix where entries are "true" above diagonal
          %  and "false" elsewhere, to allow extraction of 
          %  off-diagonal elements
  dimaway = dimaway(flag) ;

end ;
ndimaway = length(dimaway) ;



%  use estimated d0 if needed
%
if d0 == 0 ;
  sigsample = std(dimaway) ;
  d0 = (d / (2 * sigsample))^2 ;
  if sigsample > d * sqrt(1/12) ;
          %  then sample is more spread, than Unif(0,d)
          %  flag this by setting d0 to 1 
    d0 = 1 ;
  elseif d0 > d ;
          %  then distribution too tight, use upper bound
    d0 = d ;
  end ;
end ;




%  Set up theoretical values
%
mutheory = d / 2 ;
sigtheory = d / (2 * sqrt(d0)) ;
lefttheory = phi((0 - mutheory) / sigtheory) ;
righttheory = phi((d - mutheory) / sigtheory) ;
rangetheory = righttheory - lefttheory ;




%  Plot Direct Dimensions Away (if needed)
%
if  viplot(1) == 1  ;

  if iscreenwrite == 1 ;
    disp('    ViewDimAway working on direct graphic') ;
  end ;


  if nplot > 1 ;    %  then doing multiple graphics
    fighand = gcf ;
    if ~isempty(vfigpos) ;
      set(fighand,'Position',vfigpos)
    end ;
    diraxh = subplot(nplot,1,1) ;
  end ;


  dirleft = 0 ;
  dirright = d ;

  [dirkde,dirxgrid] = gpkde(dimaway,0,[dirleft; dirright]) ;
          %  0 for data driven h

  dirtheory = (1 / (sqrt(2 * pi) * sigtheory)) .* ...
                     exp(-((dirxgrid - mutheory) / sigtheory).^2 / 2) ;
          %  Unadjusted Gaussian density
  dirtheory = dirtheory / rangetheory ;
          %  Gaussian density adjusted for truncation

  dirbottom = 0 ;
  dirtop = 1.05 * max([dirkde; dirtheory]) ;


  %  compute jitter plot
  %
  if ndataoverlay == 1 ;
    ndo = min(ndimaway,1000) ;
  elseif ndataoverlay == 2 ;
    ndo = ndimaway ;
  else ;
    ndo = min(ndimaway,ndataoverlay) ;
  end ;

  if ndo < ndimaway ;    %  then need to subsample
    [temp,randperm] = sort(rand(ndimaway,1)) ;
            %  randperm is a random permutation of 1,2,...,ndimaway
    vkeep = randperm(1:ndo) ;
            %  indices of elements to keep for display
    mdatol = dimaway(vkeep,:) ;
  else ;    %  overlay full data set
    mdatol = dimaway ;
  end ;

  rhts = dirbottom + (.7 + .2 * rand(ndo,1)) * (dirtop - dirbottom) ;
              %  random heights


  dirh = plot(mdatol,rhts,'m.', ...
              dirxgrid,dirtheory,'g-', ...
              dirxgrid,dirkde,'m-') ;
    axis([dirleft,dirright,dirbottom,dirtop]) ;
    set(dirh,'LineWidth',2) ;
    if isempty(titlestr) ;
      dirtitstr = ['Dimensions Away'] ;
    else ;
      dirtitstr = [titlestr ', Dim. Away'] ;
    end ;
    dirtith = title(dirtitstr) ;



end ;    %  of Direct Plot IF block






%  Plot Probability Integral Transform (if needed)
%
if  viplot(2) == 1  ;

  if iscreenwrite == 1 ;
    disp('    ViewDimAway working on P. I. T. graphic') ;
  end ;


  if nplot > 1 ;    %  then doing multiple graphics
    diraxh = subplot(nplot,1,2) ;
  end ;


  pitdimaway = phi((dimaway - mutheory) / sigtheory) ;
          %  Probability Integral Transformed dimensions away
  pitdimaway = (pitdimaway - lefttheory) / rangetheory ;
          %  PIT version, adjusted for truncation


  pitleft = 0 ;
  pitright = 1 ;
  pitrange = pitright - pitleft ;

  padpitdimaway = [-pitdimaway; ...
                   pitdimaway; ...
                   (2 - pitdimaway) ] ;
  padpitleft = pitleft - pitrange ;
  padpitright = pitright + pitrange ;
          %  padded version, to avoid boundary effects

  hpit = bwsjpib(padpitdimaway,[padpitleft; padpitright; 601]) ;
          %  get data driven bandwidth

  [padpitkde,padpitxgrid] = gpkde(padpitdimaway,hpit, ...
                     [padpitleft; padpitright; 601]) ;
  pitkde = 3 * padpitkde(201:401) ;
  pitxgrid = padpitxgrid(201:401) ;
          %  cut off stuff for boundary adjustment

  pittheory = ones(201,1) ;

  pitbottom = 0 ;
  pittop = 1.05 * max([pitkde; pittheory]) ;


  %  compute jitter
  %
  if nplot == 1 ;
    if ndataoverlay == 1 ;
      ndo = min(ndimaway,1000) ;
    elseif ndataoverlay == 2 ;
      ndo = ndimaway ;
    else ;
      ndo = min(ndimaway,ndataoverlay) ;
    end ;
  end ;

  if ndo < ndimaway ;    %  then need to subsample
    [temp,randperm] = sort(rand(ndimaway,1)) ;
            %  randperm is a random permutation of 1,2,...,ndimaway
    pitvkeep = randperm(1:ndo) ;
            %  indices of elements to keep for display
    mdatol = pitdimaway(pitvkeep,:) ;
  else ;    %  overlay full data set
    mdatol = pitdimaway ;
  end ;

  if pittop <= 2 ;    %  then put jitter at bottom
    pitjitfrac = .1 ;
  else ;    %  then put jitter at top
    pitjitfrac = .7 ;
  end ;

  rhts = pitbottom + (pitjitfrac + .2 * rand(ndo,1)) * ...
                                   (pittop - pitbottom) ;
              %  random heights

  pith = plot(mdatol,rhts,'m.', ...
              pitxgrid,pittheory,'g-', ...
              pitxgrid,pitkde,'m-') ;
    axis([pitleft,pitright,pitbottom,pittop]) ;
    set(pith,'LineWidth',2) ;
    if  nplot > 1  |  isempty(titlestr)  ;
      pittitstr = ['Dim. Away, Subspace Gaussian Prob. Int. Trans.' ...
                          ' dim = ' num2str(d0)] ;
    else ;
      pittitstr = [titlestr ', ' num2str(d0) '-d Gauss PIT'] ;
    end ;
    pittith = title(pittitstr) ;



end ;    %  of P. I. T. Plot IF block






%  Plot Q-Q version (if needed)
%
if  viplot(3) == 1  ;

  if iscreenwrite == 1 ;
    disp('    ViewDimAway working on Q-Q plot') ;
  end ;


  if nplot > 1 ;    %  then doing multiple graphics
    diraxh = subplot(nplot,1,3) ;
  end ;


  qda = sort(dimaway) ;
          %  put in increasing order

  qgaussgrid = (1:ndimaway)' / (ndimaway + 1) ;
  qgauss = phiinv(lefttheory + rangetheory * qgaussgrid) ;
  qgauss = mutheory + qgauss * sigtheory ;



  qqleft = min([qda(1); qgauss(1)]) ;
  qqright = max([qda(ndimaway); qgauss(ndimaway)]) ;
  qqbottom = qqleft ;
  qqtop = qqright ;


  qqh = plot([qqleft;qqright],[qqbottom;qqtop],'g-', ...
              qda,qgauss,'m-') ;
    axis([qqleft,qqright,qqbottom,qqtop]) ;
    axis square ;
    set(qqh,'LineWidth',2) ;
    if  nplot > 1  |  isempty(titlestr)  ;
      qqtitstr = ['Subspace Gaussian Q-Q plot, dim = ' num2str(d0)] ;
    else ;
      qqtitstr = [titlestr ', ' num2str(d0) '-d Gauss Q-Q'] ;
    end ;
    qqtith = title(qqtitstr) ;


end ;    %  of QQ Plot IF block







%  Do main movie construction
%
if imovie == 1 ;

  if iscreenwrite == 1 ;
    disp('    ViewDimAway Working on Movie') ;
  end ;


  vdi = (dstart:dend)' ;


  mmovie(:,1) = getframe(fighand) ;



  %  Loop through and change
  %
  nframe = length(vdi) ;
  for iframe = 2:nframe ;


    di = vdi(iframe) ;

    mutheory = d / 2 ;
    sigtheory = d / (2 * sqrt(di)) ;
    lefttheory = phi((0 - mutheory) / sigtheory) ;
    righttheory = phi((d - mutheory) / sigtheory) ;
    rangetheory = righttheory - lefttheory ;


    if viplot(1) == 1 ;

      dirtheory = (1 / (sqrt(2 * pi) * sigtheory)) .* ...
                     exp(-((dirxgrid - mutheory) / sigtheory).^2 / 2) ;
          %  Unadjusted Gaussian density
      dirtheory = dirtheory / rangetheory ;
          %  Gaussian density adjusted for truncation

      set(dirh(2),'YData',dirtheory) ;
          %  reset data from green theoretical curve

    end ;


    if viplot(2) == 1 ;

      pitdimaway = phi((dimaway - mutheory) / sigtheory) ;
          %  Probability Integral Transformed dimensions away
      pitdimaway = (pitdimaway - lefttheory) / rangetheory ;
          %  PIT version, adjusted for truncation

      padpitdimaway = [-pitdimaway; ...
                       pitdimaway; ...
                       (2 - pitdimaway) ] ;
      padpitleft = pitleft - pitrange ;
      padpitright = pitright + pitrange ;
          %  padded version, to avoid boundary effects

      hpit = bwsjpib(padpitdimaway,[padpitleft; padpitright; 601]) ;
          %  get data driven bandwidth

      [padpitkde,padpitxgrid] = gpkde(padpitdimaway,hpit, ...
                     [padpitleft; padpitright; 601]) ;
      pitkde = 3 * padpitkde(201:401) ;
          %  cut off stuff for boundary adjustment

      set(pith(3),'YData',pitkde) ;
          %  reset data from magenta empirical curve


      if ndo < ndimaway ;    %  then need to subsample
        mdatol = pitdimaway(pitvkeep,:) ;
      else ;    %  overlay full data set
        mdatol = pitdimaway ;
      end ;

      set(pith(1),'XData',mdatol) ;
          %  reset data from magenta jittered data


      pittitstr = ['Dim. Away, Subspace Gaussian Prob. Int. Trans.' ...
                          ' dim = ' num2str(di)] ;
      set(pittith,'String',pittitstr)


    end ;


    if viplot(3) == 1 ;

      qgauss = phiinv(lefttheory + rangetheory * qgaussgrid) ;
      qgauss = mutheory + qgauss * sigtheory ;

      set(qqh(2),'YData',qgauss) ;
          %  reset data from magenta empirical curve

      qqtitstr = ['Subspace Gaussian Q-Q plot, dim = ' num2str(di)] ;
      set(qqtith,'String',qqtitstr)


    end ;



    mmovie(:,iframe) = getframe(fighand) ;

  end ;    %  of iframe loop, to make movie



  %  Reorder frames, and replay
  %
  vmorder = [(1:nframe),((nframe-1):-1:2)] ;
  vmorder = vec2mat(vmorder,nrif) ;
  vmorder = reshape(vmorder,1,size(vmorder,2)*nrif) ;
          %  replicate frame number, NRIF times

  mmovie = mmovie(:,vmorder') ;
          %  reorder frames, to loop back to beginning

  movie(fighand,mmovie,nrepeat) ;
          %  Play movie on screen


end ;





%  Save results (if needed)
%
if ~isempty(savestr) ;     %  then save results

  if iscreenwrite == 1 ;
    disp('    ViewDimAway saving results') ;
  end ;


  if imovie == 0 ;     %  then save as postscript file

    if  nplot == 1  ;
      orient landscape ;
    else ;
      orient tall ;
    end ;

    eval(['print -dpsc ' savestr '.ps']) ;

    writestr = 'postscript file' ;

  elseif imovie == 1 ;    %  then save as mpeg file

    comap = jet ;
          %  matlab's current default colormap

    mpgwrite(mmovie,comap,savestr,movopt) ;

    writestr = 'movie' ;

  else ;

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from ViewDimAway.m:   !!!') ;
    disp('!!!   Bad imovie parameter,       !!!') ;
    disp('!!!   Terminating Execution       !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    return ;

  end ;


  if iscreenwrite == 1 ;
    disp(['    ViewDimAway finished save of ' writestr]) ;
    disp('  ') ;
  end ;

end ;



