function [Yfit, Blocks] = monreggpav(X,Y,W);
% MONREGGPAV Monotonic regression using GPAV algorithm.
%     [Yfit, Blocks] = MONREGGPAV(X,Y,W)
%     [Yfit, Blocks] = MONREGGPAV(X,Y)
% 
%     Input:
%     * X  matrix of observations presented 
%           by one observation in each row, size [N x p]
%     * Y  response vector, size [N x 1]
%     * W  vector of weights, size [N x 1]. By default W=ones(N,1)
%     Here N is a number of observations, p - number of explanatory
%     variables
% 
%     Output:
%     * Yfit fitted response which is nondecreasing 
%           in all explanatory variables, size [N x 1]
%     * Blocks cell array containing the partitioning of indices 1..N into
%     blocks. Some cells are [0] which should be ignored.
% 
%     Last modified: 2008-09-04
%     Authors: Burdakov O., Grimvall A., Sysoev O.
%     Reference: O. Burdakov, O. Sysoev, A. Grimvall and M. Hussian. 
%                "An O(n2) algorithm for isotonic regression."
%                 In: G. Di Pillo and M. Roma (Eds) Large-Scale Nonlinear Optimization.
%                 Series: Nonconvex Optimization and Its Applications,
%                 Springer-Verlag, (2006) 83, pp. 25-33


if (~exist('W'))
    W=ones(size(X,1),1);
end
data=[X,Y,W];

[dimR, dimC] = size(data);

N=dimR;
p=dimC-2;

[data,order] = sortrows(data);
sortX = data(:,1:p);
vA = data(:, p+1);
xW = data(:, p+2);


adj=logical(sparse(N,N));

sortX=sortX';

degree=uint16(zeros(N,1));

for i=2:N
    pred = false(1,N);
    pred(1:i-1)=min(sortX(:,1:(i-1))<=repmat(sortX(:,i),1,i-1));
    degree(i)=sum(pred);
    p=find(pred);
    [pp,temp]=find(adj(:,p));
    pred(pp)=0;
    adj(:,i)=pred';
end;


[temp,degree]=sort(degree);
adj=adj(degree,degree);
  

Yfit=zeros(N,1);

vA=vA(degree);
xW=xW(degree);
[YfitN, vAN, BlocksN] = GPAV(adj,vA,xW);

degreem(degree)=1:N;
ordm(order)=1:N;
peresm=degreem(ordm);

for i=1:size(BlocksN,2)
    if (BlocksN{i}~=0)
        BlocksN{i}=degree(order(BlocksN{i}))';
    end;
end;
Yfit=YfitN(peresm);
[Blocks{1:N}]=deal(BlocksN{peresm});


