clear; clc; close all;
rng(0);
utility = struct;

% agent's utility (supposedly unknown)
alpha = 3;
u = @(x) (1-exp(-alpha*x))/(1-exp(-alpha));
invu = @(y) -log(1-y*(1-exp(-alpha)))/alpha;

f = figure; hold on;
noise_all = [1/6,1/3,1/2];
lottery_all = {};
cnt = 0;

% preparing the lotteries
n = 20;
for j = 1:length(noise_all)
noise = noise_all(j);
lottery = struct;
for i = 1:n
    lottery(i).payoffs = [(i-1)/(n-1),1];
    lottery(i).probabilities = [1,1];
    lottery(i).probabilities = ...
        lottery(i).probabilities/(sum(lottery(i).probabilities));    
    lottery(i).certaintyequivalent = ...
        invu(u(lottery(i).payoffs)*lottery(i).probabilities');
    lottery(i).certaintyequivalent = ...
        lottery(i).certaintyequivalent*(1+sin(pi*i/3)*noise);
    lottery(i).certaintyequivalent = ...
        min([lottery(i).certaintyequivalent,max(lottery(i).payoffs)]);
    lottery(i).certaintyequivalent = ...
        max([lottery(i).certaintyequivalent,min(lottery(i).payoffs)]);
end
lottery_all{j} = lottery; %#ok<SAGROW>

% non-parameteric elicitation ---------------------------------------------
[x,y,obj_np] = uhat(lottery); 

% computing coefficient of distributional robustness
codr = (y(2)-y(1))/(x(2)-x(1)); 

cnt = cnt + 1;
utility(cnt).x = x;
utility(cnt).y = y;
utility(cnt).sse = obj_np;
utility(cnt).codr = codr;

% parametric elicitation --------------------------------------------------
% assume here that u(x) = p - qexp(-rx) with
% 0 = u(0) and 1 = u(1) 
beta_range = 1e-3:1e-3:2*alpha;
obj = zeros(size(beta_range));
for s = 1:length(beta_range)
    beta = beta_range(s);
    for i = 1:n
        u2 = @(x) (1-exp(-beta*x))/(1-exp(-beta));
        obj(s) = obj(s) + ...
            (lottery(i).probabilities*u2(lottery(i).payoffs)' - ...
            u2(lottery(i).certaintyequivalent))^2;
    end
end
[~,idx] = min(obj);
beta = beta_range(idx);
obj_p = obj(idx);

fprintf('>> alpha                         = %.5f\n', alpha);
fprintf('>> perturbation magnitude        = %.5f\n', noise);
fprintf('>> parametric sum-of-square      = %.5f\n', obj_p);
fprintf('>> non-parametric sum-of-square* = %.5f\n', obj_np);
fprintf('>>                coeff robust   = %.5f\n', codr);

[noise_numerator, noise_denominator] = rat(noise);
plot(x, y, '-*', 'linewidth', 2, 'Displayname', ...
    ['\eta = ', num2str(noise_numerator), '/', num2str(noise_denominator)]);
grid on;
box on;
xlabel('$$x$$', 'Interpreter','latex');
ylabel('$$\hat{u}(x)$$','Interpreter','latex');

end

lgd = legend('show');
lgd.Location = 'southeast';
title("$\alpha = $ " + alpha, 'Interpreter', 'Latex');
plot(0:0.01:1, u(0:0.01:1), '-.', 'Color', [1 1 1]/2, 'LineWidth', 1.5, 'DisplayName', '\it u(x)');
hold off;

save(sprintf('EXP2_ALPHA=%d.mat',alpha), 'utility', 'noise_all', ...
    'lottery_all', 'u', 'invu');

