function [yhat,vaf_yhat] = monotonic_regression_dykstra(x,y)

% MONOTONIC_REGRESSION_DYKSTRA returns a vector of values in YHAT that are
% weakly monotonic with respect to the values in an independent input
% vector X, and minimize the least-squares loss to the values in the
% dependent input vector Y.  The variance-accounted-for in Y from YHAT is
% given by VAF_YHAT.
%
% syntax: [yhat,vaf_yhat] = monotonic_regression_dykstra(x,y)
%
% The least-squares optimization is carried out through Dykstra's method of
% iterative projection.

n = length(x);

[xyord,ord] = sortrows([x(:),y(:)]);

indexord(ord) = 1:n;

workspace = zeros(n-1,1);

yhatfit = xyord(:,2);

crit = 1.0;

while (crit >= 1.0e-006)
    
    crit = 0.0;
    index = 0;
    
    for jone = 1:(n-1)
        
        p1 = yhatfit(jone);
        p2 = yhatfit(jone+1);
        
        yhatfit(jone) = yhatfit(jone) + workspace(index+1);
        yhatfit(jone+1) = yhatfit(jone+1) - workspace(index+1);
        
        if((xyord(jone,1) < xyord(jone+1,1)) && ...
                (yhatfit(jone) <= yhatfit(jone+1)))
            
            workspace(index+1) = 0.0;
        end
        
       if((xyord(jone,1) < xyord(jone+1,1)) && ...
                (yhatfit(jone) > yhatfit(jone+1)))
            
            ave = (yhatfit(jone) + yhatfit(jone+1))/2;
            diff = yhatfit(jone) - ave;
            
            workspace(index+1) = diff;
            
            yhatfit(jone) = ave;
            yhatfit(jone+1) = ave;
       end
       
       index = index + 1;
       
       crit = crit + abs(p1 - yhatfit(jone)) + abs(p2 - yhatfit(jone+1));
       
    end
    
end
    
    aver = (sum(xyord(:,2)))/n;
    
    numer = sum((xyord(:,2) - yhatfit(:)).^2);
    denom = sum((xyord(:,2) - aver).^2);
    
    vaf_yhat = 1 - (numer/denom);
    
    yhat = yhatfit(indexord);
    
    plot(xyord(:,1),xyord(:,2),'k.')
    hold on
    plot(xyord(:,1),yhatfit(:),'r*')
    hold on
    plot(xyord(:,1),yhatfit(:),'g-')
    
    
            