%
% Takashi Goda's model testcase from our paper
% 'Decision-making under uncertainty: using
%  MLMC for efficient estimation of EVPPI'
% Statistics and Computing, 29(4):739-751, 2019.
%

function nested

close all; clear all;

%
% MLMC treatment
%

addpath('..');

N     = 200000; % samples for convergence tests
L     = 10;     % levels for convergence tests 

N0    = 1000;   % initial samples on coarse levels
Lmin  = 2;      % minimum refinement level
Lmax  = 20;     % maximum refinement level

Eps    = [ 0.0001 0.0002 0.0005 0.001 0.002 ];

filename = 'nested';
fp       = fopen([filename '.txt'],'w');
mlmc_test(@test3_l, N,L, N0,Eps,Lmin,Lmax, fp);
fclose(fp);

%
% plot results
%

nvert = 3;
mlmc_plot(filename, nvert);

if(nvert==1)
  figure(1)
  print('-deps2',[filename 'a.eps'])
  figure(2)
  print('-deps2',[filename 'b.eps'])
else
  print('-deps2',[filename '.eps'])
end

%-------------------------------------------------------
%
% level l estimator
%

function [sums cost] = test3_l(l,N)

% global option

nf = 2^l;
nc = nf/2;

sums(1:6) = 0;

for N1 = 1:1000:N
  N2 = min(1000,N-N1+1);

% level 0

  if (l==0)
    Pf = zeros(1,N2);
    dP = Pf;

  else

% level l>0, with antithetic sampler

    X  = randn(1,N2);  % outer variables
    X  = repmat(X,nf,1);
    Y  = randn(nf,N2); % inner variables
    
    f2 = X + Y;

    ave = sum(max(0,f2),1)/nf;
    Pf  = ave - max(0,sum(f2,1)/nf);
    Pc  = ave - 0.5*( max(0,sum(f2(   1:nc,:),1)/nc) ...
                    + max(0,sum(f2(nc+1:nf,:),1)/nc) );
    dP = Pf - Pc;
  end

  sums(1) = sums(1) + sum(dP);
  sums(2) = sums(2) + sum(dP.^2);
  sums(3) = sums(3) + sum(dP.^3);
  sums(4) = sums(4) + sum(dP.^4);
  sums(5) = sums(5) + sum(Pf);
  sums(6) = sums(6) + sum(Pf.^2);
end

cost = N*nf;
