% Example 3.3
% This illustrates the concepts and procedure for isotonic regression.

% Enter the original data.
dissim = [2.1 3 2.4 1.7 3.9 3.2];
dists = [1.6 4.5 5.7 3.3 4.3 1.3];
N = length(dissim);
% Now re-order the dissimilarities.
[dissim,ind] = sort(dissim);
% Now impose the same order on the distances.
dists = dists(ind);
% Now find the cumulative sums of the distances.
D = cumsum(dists);
% Add the origin as the first point.
D = [0 D];

% Now find the slope of these.
slope = D(2:end)./(1:N);
% Find the points on the convex minorant by looking
% for smallest slopes.
i = 1;
k = 1;
while i <= N
    val = min(slope(i:N));
    minpt(k) = find(slope == val);
    i = minpt(k) + 1;
    k = k + 1;
end

% Now that we have all of the minorant points that divide into blocks, 
% the disparities are the averages of the distances over those blocks.
j = 1;
for i = 1:length(minpt)
    dispars(j:minpt(i)) = mean(dists(j:minpt(i)));
    j = minpt(i) + 1;
end

% Now for the plot.
cvm = [0 minpt];
plot(0:6,D, '-ko',cvm,D(cvm + 1),'k--')
xlabel('i'),ylabel('D_i')
legend({'Cumulative Sum','Greatest Convex Minorant'})


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Better way to find points on the convex hull, written by Tom Lane, THe
% MathWorks.

% Set up with fake data
rand('state',0);
dists = [1 0 1 .9]; N = length(dists);
dists = abs(rand(1,10) + (1:10)/15); N = length(dists);
D = cumsum(dists);
D = [0 D];
clear minpt newminpt

% Code as in text
i = 1;
k = 1;
slope = D(2:end) ./ (1:N);
while (i<=N)
    val = min(slope(i:N));
    minpt(k) = find(slope == val);
    i = minpt(k) + 1;
    k = k+1;
end
minpt

% Tom's correction
k = 1;
offset = 0;
todo = 1:N;
while ~isempty(todo)
    slope = (D(1+todo)-D(1+offset))./(todo-offset); % slope from prev point
    val = min(slope);                % find min slope of remaining points
    nextpt = min(find(slope == val));% get next point
    newminpt(k) = offset + nextpt;   % get its distance from zero
    offset = newminpt(k);            % compute new slopes from this origin
    todo = offset+1:N;               % remaining points to consider
    k = k+1;
end
newminpt
plot(0:N,D,'gx',[0 minpt],D(1+[0 minpt]),'bo',[0 newminpt],D(1+[0 newminpt]),'rs:')
legend('All points','Hull using text method','Hull using my method',2)
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Alternative way to find the disparaties by Tom Lane, The MathWorks, Inc.

%% Page 80 - alternative way to compute disparities
minpt = newminpt;

% As in text, looping and taking means of groups of points
j = 1;
for i = 1:length(minpt)
    dispars(j:minpt(i)) = mean(dists(j:minpt(i)));
    j = minpt(i) + 1;
end
dispars

% By interpolating along convex boundary and then differencing
cumdispars = interp1([0,minpt], [0,D(1+minpt)], 0:N);
dispars = diff(cumdispars)
