function h = qqplot2(x,y,pvec)
% Bahman Javadi
% I have changed the standard function 'qqplot' to plot in a rectangle with x=y
% line
% INRIA 2008-2009
% 
%QQPLOT Display an empirical quantile-quantile plot.
%   QQPLOT(X) makes an empirical QQ-plot of the quantiles of
%   the data set X versus the quantiles of a standard Normal distribution.
%
%   QQPLOT(X,Y) makes an empirical QQ-plot of the quantiles of
%   the data set X versus the quantiles of the data set Y.
%
%   H = QQPLOT(X,Y,PVEC) allows you to specify the plotted quantiles in 
%   the vector PVEC. H is a handle to the plotted lines. 
%
%   When both X and Y are input, the default quantiles are those of the 
%   smaller data set.
%
%   The purpose of the quantile-quantile plot is to determine whether
%   the sample in X is drawn from a Normal (i.e., Gaussian) distribution,
%   or whether the samples in X and Y come from the same distribution
%   type.  If the samples do come from the same distribution (same shape),
%   even if one distribution is shifted and re-scaled from the other
%   (different location and scale parameters), the plot will be linear.
%
%   Use the data cursor to read precise observation values and where 
%   they project to on the reference line.  The variable's observation 
%   numbers are displayed when they have have no missing values and are no 
%   longer than the other variable.
%
%   See also NORMPLOT.

%   Copyright 1993-2008 The MathWorks, Inc. 
%   $Revision: 2.13.2.2 $  $Date: 2008/02/29 13:12:18 $

if nargin == 1
    if size(x,1)==1
        x = x';
    end
    [y,origindy]  =  sort(x);
    x = plotpos(y);
    origindx = [];
    x = norminv(x);
    xx = x;
    yy = y;
else
    if size(x,1)==1
        x = x';
    end
    if size(y,1)==1
        y = y';
    end
    n = -1;
    if nargin < 3
        nx = sum(~isnan(x));
        if (length(nx) > 1)
            nx = max(nx);
        end
        ny = sum(~isnan(y));
        if (length(ny) > 1)
            ny = max(ny);
        end
        n = min(nx, ny);
        pvec = 100*((1:n) - 0.5) ./ n;
    end
    if size(x,1)==n
        xx = zeros(size(x));
        origindx = zeros(size(x));
        nancols = find(any(isnan(x),1));
        fullcols = find(all(~isnan(x),1));
        [xx(:,fullcols),origindx(:,fullcols)] = sort(x(:,fullcols));
        xx(:,nancols) = prctile(x(:,nancols),pvec);
        if size(x,2)==1 && size(y,2)~=1
            origindx = repmat(origindx,1,size(y,2));
        end

    else
        xx = prctile(x,pvec);
        origindx = [];
    end
    if size(y,1)==n
        yy = zeros(size(y));
        origindy = zeros(size(y));
        nancols = find(any(isnan(y),1));
        fullcols = find(all(~isnan(y),1));
        [yy(:,fullcols),origindy(:,fullcols)] = sort(y(:,fullcols));
        yy(:,nancols) = prctile(y(:,nancols),pvec);
        if size(y,2)==1 && size(x,2)~=1
            origindy = repmat(origindy,1,size(x,2));
        end
    else
        yy = prctile(y,pvec);
        origindy = [];
    end
end

q1x = prctile(x,25);
q3x = prctile(x,75);
q1y = prctile(y,25);
q3y = prctile(y,75);
qx = [q1x; q3x];
qy = [q1y; q3y];


dx = q3x - q1x;
dy = q3y - q1y;
slope = dy./dx;
centerx = (q1x + q3x)/2;
centery = (q1y + q3y)/2;
maxx = max(x);
minx = min(x);
maxy = centery + slope.*(maxx - centerx);
miny = centery - slope.*(centerx - minx);

mx = [minx; maxx];
my = [miny; maxy];

newplot();
%hrefends = line(mx,my,'LineStyle','-.','Marker','none');
%hrefmid = line(qx,qy,'LineStyle','-','Marker','none');
%added by Bahman
qxx = [0:1:numel(xx)];

hrefends = line(qxx,qxx,'LineStyle','-.','Marker','none');
hrefmid = line(qxx,qxx,'LineStyle','-','Marker','none');

hdat = line(xx,yy,'LineStyle','none','Marker','+');

%{
%---------------------- log scale
% added by Bahman
hrefends = loglog(qxx,qxx,'LineStyle','-.','Color','r');
hold on
hdat = loglog(xx,yy,'LineStyle','none','Marker','+');
return;
%-------------------------------
%}

if length(hdat)==1
    set(hdat,'MarkerEdgeColor','b');
    set([hrefends,hrefmid],'Color','r');
end
if nargout>0
    h = [hdat;hrefmid;hrefends];
end

for i=1:length(hdat)
    % Set custom data cursor on data
    hB = hggetbehavior(hdat(i),'datacursor');
    set(hB,'UpdateFcn',@qqplotDatatipCallback);
    % Disable datacursor on reference lines
    hB = hggetbehavior(hrefends(i),'datacursor');
    set(hB,'Enable',false);
    hB = hggetbehavior(hrefmid(i),'datacursor');
    set(hB,'Enable',false);
    if length(hdat)>1
        setappdata(hdat(i),'group',i);
    end
    if ~isempty(origindx)
        setappdata(hdat(i),'origindx',origindx(:,i));
    end
    if ~isempty(origindy)
        setappdata(hdat(i),'origindy',origindy(:,i));
    end
end

if nargin == 1
   xlabel('Standard Normal Quantiles')
   ylabel('Quantiles of Input Sample')
   title ('QQ Plot of Sample Data versus Standard Normal')
else
   xlabel('X Quantiles');
   ylabel('Y Quantiles');
end

%===================== helper function ====================
function [pp,n] = plotpos(sx)
%PLOTPOS Compute plotting positions for a probability plot
%   PP = PLOTPOS(SX) compute the plotting positions for a probability
%   plot of the columns of SX (or for SX itself if it is a vector).
%   SX must be sorted before being passed into PLOTPOS.  The ith
%   value of SX has plotting position (i-0.5)/n, where n is
%   the number of rows of SX.  NaN values are removed before
%   computing the plotting positions.
%
%   [PP,N] = PLOTPOS(SX) also returns N, the largest sample size
%   among the columns of SX.  This N can be used to set axis limits.

[n, m] = size(sx);
if n == 1
   sx = sx';
   n = m;
   m = 1;
end

nvec = sum(~isnan(sx));
pp = repmat((1:n)', 1, m);
pp = (pp-.5) ./ repmat(nvec, n, 1);
pp(isnan(sx)) = NaN;

if (nargout > 1)
   n = max(nvec);  % sample size for setting axis limits
end

%===================== callback ====================
function datatipTxt = qqplotDatatipCallback(obj,evt)

target = get(evt,'Target');
pos = get(evt,'Position');
ind = get(evt,'DataIndex');

group = getappdata(target,'group');
origindx = getappdata(target,'origindx');
origindy = getappdata(target,'origindy');


datatipTxt = {
    ['x: ',num2str(pos(1))],...
    ['y: ',num2str(pos(2))]...
    };

if ~isempty(group) || ~isempty(origindx) || ~isempty(origindy)
    datatipTxt{end+1} = '';
end

if ~isempty(group)
    datatipTxt{end+1} = ['Group: ',num2str(group)];
end
if ~isempty(origindx)
    dat = origindx(ind);
    if dat~=0
        datatipTxt{end+1} = ['X observation: ',num2str(dat)];
    end
end
if ~isempty(origindy)
    dat = origindy(ind);
    if dat~=0
        datatipTxt{end+1} = ['Y observation: ',num2str(dat)];
    end
end
