clear; close all;

f = figure; hold on;
color = get(gca,'colororder');

for alpha = 5:-1:3

% utility elicitation 
load(sprintf('EXP1_ALPHA=%d.mat',alpha));
x = utility(2).x;
y = utility(2).y;   
lottery = lottery_all{2};

% compute sum of squared errors before reduction
ucal = @(in) ...
    min(((y(2:end)-y(1:end-1))*in + ...
          x(2:end).*y(1:end-1) - ...
          x(1:end-1).*y(2:end))./(x(2:end)-x(1:end-1)));
sse = 0;
for l = 1:length(lottery)
    e = 0;
    for ll = 1:length(lottery(l).payoffs)
        e = e + lottery(l).probabilities(ll)*ucal(lottery(l).payoffs(ll));
    end
    e = e - ucal(lottery(l).certaintyequivalent);
    sse = sse + e^2;
end
fprintf('>> alpha            = %d\n', alpha)
fprintf('>> Original SSE     = %.6f\n', sse); % should = utility(2).sse

SIZE = 4;
while length(x) > SIZE
    c = zeros(1,length(x));
    c(1) = -inf;
    c(end) = -inf;
    for i = 2:length(x)-1
        c(i) = (y(i+1)-y(i))/(x(i+1)-x(i));
        c(i) = c(i)-(y(i)-y(i-1))/(x(i)-x(i-1));
    end
    [~,idx] = max(c);
    x = x([1:idx-1,idx+1:end]);
    y = y([1:idx-1,idx+1:end]);
end

figure(f);
ax = gca(f);
ax.FontSize = 14;

% before reduction
p = plot(utility(2).x, utility(2).y, '-', 'Color', color(alpha-2,:), ...
    'Linewidth', 1.5, 'HandleVisibility', 'off');

% after reduction
plot(x, y, '-', 'Color', color(alpha-2,:), 'Linewidth', 1.5, ...
    'Displayname', ['\alpha = ', num2str(alpha)]);

% reduction gap
fill([utility(2).x, x(end:-1:1)], [utility(2).y, y(end:-1:1)], ...
    color(alpha-2,:), 'HandleVisibility', 'off');
grid on;
box on;
axis([0 1 0 1]);
xlabel('$$x$$','Interpreter','latex','FontSize',22);
ylabel('$$\hat{u}(x)$$','Interpreter','latex','FontSize',22);
titletext = sprintf("|\\mathcal{K}| = %d", SIZE); 
title(strcat("$", titletext, "$"),'Interpreter','latex','FontSize',22);

% compute the sum of squared errors after reduction
ucal = @(in) ...
    min(((y(2:end)-y(1:end-1))*in + ...
          x(2:end).*y(1:end-1) - ...
          x(1:end-1).*y(2:end))./(x(2:end)-x(1:end-1)));
sse = 0;
for l = 1:length(lottery)
    e = 0;
    for ll = 1:length(lottery(l).payoffs)
        e = e + lottery(l).probabilities(ll)*ucal(lottery(l).payoffs(ll));
    end
    e = e - ucal(lottery(l).certaintyequivalent);
    sse = sse + e^2;
end
fprintf('>> New SSE (approx) = %.6f\n', sse);

% miqp-based approach
[x,y] = uhat2(lottery,SIZE);
ucal = @(in) ...
    min(((y(2:end)-y(1:end-1))*in + ...
          x(2:end).*y(1:end-1) - ...
          x(1:end-1).*y(2:end))./(x(2:end)-x(1:end-1)));
sse = 0;
for l = 1:length(lottery)
    e = 0;
    for ll = 1:length(lottery(l).payoffs)
        e = e + lottery(l).probabilities(ll)*ucal(lottery(l).payoffs(ll));
    end
    e = e - ucal(lottery(l).certaintyequivalent);
    sse = sse + e^2;
end
fprintf('>> New SSE (exact)  = %.6f\n', sse);

end

hold off;
lgd = legend('show');
lgd.Location = 'southeast';

