function y = weightpercent(X,p,w)
%Returns weighted percentiles of a sample given the weight vector w
% The idea is to give more emphasis in some examples of data as compared to
% others by giving more weight. For example, we could give lower weights to
% the outliers. 
% The motivation to write this function is to compute percentiles for Monte
% Carlos simulations where some simulations are very bad (in terms of goodness
% of fit between simulated and actual value) than the others and to give 
% the lower weights based on some goodness of fit criteria.
%
% USAGE:
%         y = WPRCTILE(X,p,w)
%                                             
% INPUT:
%    X -  vector or matrix of the sample data                                 
%    p -  scalar  or a vector of percent values between 0 and 100
%    w -  positive weight vector for the sample data. Length of w must be
%         equal to either number of rows or columns of X.
%         If X is matrix, same weight vector w is used for all columns (for DIM=1)
%         or for all rows (for DIM=2).
%         If the weights are equal, then  WPRCTILE is same as PRCTILE, however
%         there will be slight difference in the values of y due to different
%         implementation
%       
% OUTPUT:
%    y - percentiles of the values in X
%        When X is a vector, y is the same size as p, and y(i) contains the
%        P(i)-th percentile.
%        When X is a matrix, WPRCTILE calculates percentiles along dimension DIM
%        which is based on: if size(X,1) == length(w), DIM = 1;
%                       elseif size(X,2) == length(w), DIM = 2;
%                      
% Please note that this version of WPRCTILE will not work with NaNs values and
% planned to update in near future to handle NaNs values as missing values.
%
% EXAMPLES:
%    w = rand(1000,1);
%    y = prctile(x,[2.5 25 50 75 97.5],w);
%    % here if the size of x is 1000-by-50, then y will be size of 6-by-50
%    % if x is 50-by-1000, then y will be of the size of 50-by-6
%     

%% Input arguments check

error(nargchk(3,3,nargin))
if ~isvector(p) || numel(p) == 0
    error('wprctile:BadPercents', ...
          'P must be a scalar or a non-empty vector.');
elseif any(p < 0 | p > 100) || ~isreal(p)
    error('wprctile:BadPercents', ...
          'P must take real values between 0 and 100');
end
if ndims(X)>2
   error('wprctile:InvalidNumberofDimensions','X Must be 2D.')
end

if ~isvector(w)|| any(w<0) 
	error('wprctile:InvalidWeight', ...
          'w must vecor and values should be greater than 0');
end

% Check if there are NaN in any of the input
nans = isnan(X);
if any(nans(:)) || any(isnan(p))|| any(isnan(w))
	error('wprctile:NaNsInput',['This version of WPRCTILE will not work with ' ...
	      'NaNs values in any input and planned to update in near future to ' ...
		   'handle NaNs values as missing values.']);
end
%% Figure out which dimension WPRCTILE will work along using weight vector w

nw = length(w);
[nrows ncols] = size(X);
if nrows==nw
	dim = 1;
elseif ncols==nw
	dim = 2;
else
	error('wprctile:InvalidDimension', ...
          'length of w must be equal to either number of rows or columns of X');
end

%% Work along DIM = 1 i.e. columswise, convert back later if needed using tflag

tflag = false; % flag to note transpose
if dim==2     
   X = X';
   tflag = true;  
end
ncols = size(X,2);
np = length(p);
y = zeros(np,ncols);

% Change w to column vector and normalise it to have sum of 1
w = w(:);
w = w/sum(w);

%% Work on each column separately because of weight vector

for i=1:ncols
	temp = [X(:,i) w];
	sortedX = sortrows(temp); % sort
	cumW = cumsum(sortedX(:,2));
	xx = sortedX(:,1);
	% to avoid NaN for outside the range, the minimum or maximum values in X are
	% assigned to percentiles for percent values outside that range.
	q = [0;cumW;1];
	xxe = [xx(1); xx; xx(end)];
	y(:,i) = interp1q(q,xxe,p(:)./100);
end

%% Transpose data back for DIM = 2 to the orginal dimension of X

if tflag 
	y=y';
end