function [h,ax,BigAx,patches] = plotmatrix(varargin)
%PLOTMATRIX Scatter plot matrix.
%   PLOTMATRIX(X,Y) scatter plots the columns of X against the columns
%   of Y.  If X is P-by-M and Y is P-by-N, PLOTMATRIX will produce a
%   N-by-M matrix of axes. PLOTMATRIX(Y) is the same as PLOTMATRIX(Y,Y)
%   except that the diagonal will be replaced by HIST(Y(:,i)). 
%
%   PLOTMATRIX(...,'LineSpec') uses the given line specification in the
%   string 'LineSpec'; '.' is the default (see PLOT for possibilities).  
%
%   [H,AX,BigAx,P] = PLOTMATRIX(...) returns a matrix of handles to the
%   objects created in H, a matrix of handles to the individual subaxes
%   in AX, a handle to a big (invisible) axes which frames the
%   subaxes in BigAx, and the a matrix of handles for the histogram
%   plots in P.  BigAx is left as the CurrentAxes so that a subsequent
%   TITLE, XLABEL, or YLABEL will be centered with respect to the
%   matrix of axes.
%
%   Example:
%       x = randn(50,3); y = x*[-1 2 1;2 0 1;1 -2 3;]';
%       plotmatrix(y)

%   Clay M. Thompson 10-3-94
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.12 $  $Date: 1998/12/18 16:28:55 $

error(nargchk(1,3,nargin));
nin = nargin;

sym = '.'; % Default scatter plot symbol.
plotfun = 'plot';
dohist = 0;

if isstr(varargin{nin}),
  sym = varargin{nin};
  nin = nin - 1;
end

if nin==1, % plotmatrix(y)
  rows = size(varargin{1},2); cols = rows;
  x = varargin{1}; y = varargin{1};
  dohist = 1;
elseif nin==2, % plotmatrix(x,y)
  rows = size(varargin{2},2); cols = size(varargin{1},2);
  x = varargin{1}; y = varargin{2};
end

% Don't plot anything if either x or y is empty
if isempty(rows) | isempty(cols),
   if nargout>0, h = []; ax = []; BigAx = []; end
   return
end

if (ndims(x)==2) & any(size(x)==1), x = x(:); end
if (ndims(y)==2) & any(size(y)==1), y = y(:); end

if ndims(x)>2 | ndims(y)>2,
  error('X and Y must be 2-D.')
end
if size(x,1)~=size(y,1) | size(x,3)~=size(y,3),
  error('X and Y must have the same number of rows and pages.');
end

% Create/find BigAx and make it invisible
BigAx = newplot;
next = lower(get(BigAx,'NextPlot'));
hold_state = ishold;
set(BigAx,'Visible','off','color','none')

if any(sym=='.'),
  units = get(BigAx,'units');
  set(BigAx,'units','pixels');
  pos = get(BigAx,'Position');
  set(BigAx,'units',units);
  markersize = max(1,min(15,round(15*min(pos(3:4))/size(x,1)/max(rows,cols))));
else
  markersize = get(0,'defaultlinemarkersize');
end

% Create and plot into axes
ax = zeros(rows,cols);
pos = get(BigAx,'Position');
width = pos(3)/cols;
height = pos(4)/rows;
space = .02; % 2 percent space between axes
pos(1:2) = pos(1:2) + space*[width height];
[m,n,k] = size(y);
xlim = zeros([rows cols 2]);
ylim = zeros([rows cols 2]);
for i=rows:-1:1,
  for j=cols:-1:1,
    axPos = [pos(1)+(j-1)*width pos(2)+(rows-i)*height ...
             width*(1-space) height*(1-space)];
    findax = findobj(gcf,'Type','axes','Position',axPos);
    if isempty(findax),
      ax(i,j) = axes('Position',axPos);
      set(ax(i,j),'visible','on');
    else
      ax(i,j) = findax(1);
    end
    hh(i,j,:) = feval(plotfun,reshape(x(:,j,:),[m k]), ...
                              reshape(y(:,i,:),[m k]),sym)';
    set(hh(i,j,:),'markersize',markersize);
    set(ax(i,j),'xlimmode','auto','ylimmode','auto','xgrid','off','ygrid','off')
    xlim(i,j,:) = get(gca,'xlim');
    ylim(i,j,:) = get(gca,'ylim');
  end
end

xlimmin = min(xlim(:,:,1),[],1); xlimmax = max(xlim(:,:,2),[],1);
ylimmin = min(ylim(:,:,1),[],2); ylimmax = max(ylim(:,:,2),[],2);

% Try to be smart about axes limits and labels.  Set all the limits of a
% row or column to be the same and inset the tick marks by 10 percent.
inset = .15;
for i=1:rows,
  set(ax(i,1),'ylim',[ylimmin(i,1) ylimmax(i,1)])
  dy = diff(get(ax(i,1),'ylim'))*inset;
  set(ax(i,:),'ylim',[ylimmin(i,1)-dy ylimmax(i,1)+dy])
end
for j=1:cols,
  set(ax(1,j),'xlim',[xlimmin(1,j) xlimmax(1,j)])
  dx = diff(get(ax(1,j),'xlim'))*inset;
  set(ax(:,j),'xlim',[xlimmin(1,j)-dx xlimmax(1,j)+dx])
end

set(ax(1:rows-1,:),'xticklabel','')
set(ax(:,2:cols),'yticklabel','')
set(BigAx,'XTick',get(ax(rows,1),'xtick'),'YTick',get(ax(rows,1),'ytick'), ...
          'userdata',ax,'tag','PlotMatrixBigAx')

if dohist, % Put a histogram on the diagonal for plotmatrix(y) case
  for i=rows:-1:1,
    histax = axes('Position',get(ax(i,i),'Position'));
    [nn,xx] = hist(reshape(y(:,i,:),[m k]));
    patches(i,:) = bar(xx,nn,'hist');
    set(histax,'xtick',[],'ytick',[],'xgrid','off','ygrid','off')
  end
  patches = patches';
end

% Make BigAx the CurrentAxes
set(gcf,'CurrentAx',BigAx)
if ~hold_state,
   set(gcf,'NextPlot','replace')
end

% Also set Title and X/YLabel visibility to on and strings to empty
set([get(BigAx,'Title'); get(BigAx,'XLabel'); get(BigAx,'YLabel')], ...
 'String','','Visible','on')

if nargout~=0,
  h = hh;
end
