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

% CONVEX_CONCAVE_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, define a convex function up to the median 
% observation in X and a concave function thereafter, 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] = convex_concave_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(4*n,1);

yhatfit = xyord(:,2);

crit = 1.0;

if (mod(n,2) == 1)
    
    upper_limit = (n+3)/2;
    lower_limit = (n-1)/2;
    
end

if (mod(n,2) == 0)
    
    upper_limit = (n/2) + 1;
    lower_limit = (n/2);
    
end

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
    
    for jone = 1:(upper_limit-2)
        
        if ((xyord(jone+2,1) - xyord(jone,1)) > 1.0e-014)   
        
        p1 = yhatfit(jone);
        p2 = yhatfit(jone+1);
        p3 = yhatfit(jone+2);
        
        yhatfit(jone) = yhatfit(jone) - workspace(index+1);
        yhatfit(jone+1) = yhatfit(jone+1) + workspace(index+2);
        yhatfit(jone+2) = yhatfit(jone+2) - workspace(index+3);
        
        weight_beta = (xyord(jone+2,1) - xyord(jone+1,1))/...
            (xyord(jone+2,1) - xyord(jone,1));
        
        weight_delta = (xyord(jone+1,1) - xyord(jone,1))/...
            (xyord(jone+2,1) - xyord(jone,1));
        
        div = (weight_beta)^2 + (weight_delta)^2 + 1.0;
        
        discrep = ((-1.0)*weight_beta*yhatfit(jone) + yhatfit(jone+1) + ...
            (-1.0)*weight_delta*yhatfit(jone+2))/div;
            
       
        if (discrep <= 0.0)
            
            workspace(index+1) = 0.0; 
            workspace(index+2) = 0.0;
            workspace(index+3) = 0.0;
            
            
        end
        
        if (discrep > 0.0)
            
            workspace(index+1) = discrep*weight_beta;
            workspace(index+2) = discrep;
            workspace(index+3) = discrep*weight_delta;
            
            yhatfit(jone) = yhatfit(jone) + discrep*weight_beta;
            yhatfit(jone+1) = yhatfit(jone+1) + (-1.0)*discrep;
            yhatfit(jone+2) = yhatfit(jone+2) + discrep*weight_delta;
            
        end
        
        index = index + 3;
        
        crit = crit + abs(p1-yhatfit(jone)) + abs(p2-yhatfit(jone+1)) + ...
            abs(p3-yhatfit(jone+2));
        
        end
    end
    
    for jone = lower_limit:(n-2)
        
      if ((xyord(jone+2,1) - xyord(jone,1)) > 1.0e-014)   
        
        p1 = yhatfit(jone);
        p2 = yhatfit(jone+1);
        p3 = yhatfit(jone+2);
        
        yhatfit(jone) = yhatfit(jone) + workspace(index+1);
        yhatfit(jone+1) = yhatfit(jone+1) - workspace(index+2);
        yhatfit(jone+2) = yhatfit(jone+2) + workspace(index+3);
        
        weight_beta = (xyord(jone+2,1) - xyord(jone+1,1))/...
            (xyord(jone+2,1) - xyord(jone,1));
        
        weight_delta = (xyord(jone+1,1) - xyord(jone,1))/...
            (xyord(jone+2,1) - xyord(jone,1));
        
        div = (weight_beta)^2 + (weight_delta)^2 + 1.0;
        
        discrep = (weight_beta*yhatfit(jone) + (-1.0)*yhatfit(jone+1) + ...
            weight_delta*yhatfit(jone+2))/div;
            
       
        if (discrep <= 0.0)
            
            workspace(index+1) = 0.0; 
            workspace(index+2) = 0.0;
            workspace(index+3) = 0.0;
            
            
        end
        
        if (discrep > 0.0)
            
            workspace(index+1) = discrep*weight_beta;
            workspace(index+2) = discrep;
            workspace(index+3) = discrep*weight_delta;
            
            yhatfit(jone) = yhatfit(jone) + (-1.0)*discrep*weight_beta;
            yhatfit(jone+1) = yhatfit(jone+1) + discrep;
            yhatfit(jone+2) = yhatfit(jone+2) + (-1.0)*discrep*weight_delta;
            
        end
        
        index = index + 3;
        
        crit = crit + abs(p1-yhatfit(jone)) + abs(p2-yhatfit(jone+1)) + ...
            abs(p3-yhatfit(jone+2));
        
      end       
           
    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-')
    
    
    
            