
function [YFit, vA, Blocks] = GPAV(adj,Y,W);
%GPAV   Monotonic regression using GPAV algorithm.
%     Yfit = GPAV(adj,Y,W)
% 
%     Input:
%     * adj adjacency matrix which is sparse, upper triangular 
%       logical matrix, size [N x N]
%     * 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

N=size(adj,1);
corr=1:N;
vA=Y;
xW=W;

for q = 1:N
    Blocks{q} = [q];
end

for j=1:N
    Items = [1];
    while ~isempty(Items)
        [jMinus,temp] = find(adj(:,Blocks{j}));
        jMinus=corr(jMinus);
        Items = find(vA(jMinus) > vA(j));

        if ~isempty(Items)
            LookingIn = jMinus(Items);
            [Xmax,I] = max(vA(LookingIn));
            i = LookingIn(I(1));

            %Merging
            vA(j) = (xW(i)*vA(i) + xW(j)*vA(j))/(xW(j) + xW(i));
            xW(j) = xW(j) + xW(i);

            % Add block i into j
            Blocks{j} = [Blocks{i}, Blocks{j}];
            Blocks{i} = [0];
            corr(Blocks{j})=j;

            vA(i)=0;
        end
    end
end

YFit=vA(corr);
 





