% Generate efficiency plot.

%% DATA
log_instantaneous_acceptance = -76.5;
nsteps_list = [ 1        2        4        8       16       32       64      128      256      512     1024     2048 ];
log_mean_acceptance = [-76.4    -78.1    -80.3    -59.9    -41.1    -28.9    -15.3    -10.6     -6.8     -5.0     -3.4     -2.4];
log_relative_efficiency = [ 0.1     -2.3     -5.2     14.5     32.6     44.2     57.0     61.1     64.1     65.2     66.2     66.4 ];

log_instantaneous_acceptance = -25.7;
nsteps_list = [1        2        4        8       16       32       64      128      256      512     1024     2048];
log_mean_acceptance = [-25.7    -26.4    -27.6    -30.4    -33.7    -36.1    -18.0    -14.7     -6.1     -4.7     -3.1     -2.2];
log_mean_acceptance_error = [1.000    1.000    1.000    1.000    1.000    0.721    0.783    0.302    0.370    0.156    0.073    0.042];
log_relative_efficiency = [-0.7     -1.7     -3.5     -6.8    -10.8    -13.8      3.5      6.2     14.0     14.8     15.7     15.9 ];

% 1770 iterations
log_instantaneous_acceptance = -70.1;
log_mean_acceptance = [-70.1    -70.7    -72.9    -74.5    -60.5    -33.0    -16.9    -12.8     -8.2     -4.5     -3.1     -2.0];
log_mean_acceptance_error = [0.000    0.000    0.000    0.000    0.000    0.989    0.990    0.724    0.578    0.184    0.091    0.051];
log_relative_efficiency = [-0.6     -1.7     -4.4     -6.6      6.8     33.6     49.1     52.5     56.4     59.4     60.1     60.5];

% 4395 iterations
log_instantaneous_acceptance = -41.1;
log_mean_acceptance = [ -41.2    -42.7    -45.8    -53.2    -45.5    -24.7    -17.6    -10.1     -7.2     -4.5     -3.1     -2.1];
log_mean_acceptance_error = [1.000    1.000    1.000    1.000    1.000    1.000    0.829    0.775    0.382    0.116    0.058    0.033];
log_relative_efficiency = [ -0.7     -2.7     -6.3    -14.3     -7.2     12.9     19.3     26.1     28.4     30.4     31.1     31.4];

% 5190 iterations
log_instantaneous_acceptance = -41.3;
nsteps_list = [      1      2      4      8     16     32     64    128    256    512   1024   2048 ];
log_mean_acceptance = [ -41.3 -42.8 -46.0 -53.4 -45.7 -24.9 -17.8 -10.3  -7.1  -4.4  -3.1  -2.1 ];
log_mean_acceptance_lower = [ -81.3 -83.0 -86.6 -80.5 -61.5 -38.7 -21.1 -12.9  -8.0  -4.6  -3.2  -2.2 ];
log_mean_acceptance_upper = [ -40.2 -41.7 -44.9 -52.3 -44.6 -23.8 -16.8  -9.3  -6.5  -4.2  -3.0  -2.0 ];
log_relative_efficiency = [  -0.7  -2.7  -6.3 -14.3  -7.2  12.9  19.3  26.1  28.7  30.6  31.2  31.6 ];
log_relative_efficiency_lower = [  -0.7  -2.7  -6.3 -14.3 -23.1  -0.8  15.8  23.3  27.2  29.4  30.0  30.4 ];
log_relative_efficiency_upper = [  -0.6  -1.7  -4.3   5.7  32.0  48.8  59.3  65.7  68.7  70.7  71.3  71.6 ];

% 6085 iterations
log_instantaneous_acceptance = -41.4;
nsteps_list = [      1      2      4      8     16     32     64    128    256    512   1024   2048 ];
log_mean_acceptance = [ -41.5 -43.0 -46.2 -53.6 -45.8 -25.1 -17.8 -10.4  -7.0  -4.5  -3.1  -2.1 ];
log_mean_acceptance_lower = [ -73.4 -74.4 -77.7 -75.7 -61.7 -38.7 -20.4 -13.0  -7.8  -4.7  -3.2  -2.1 ];
log_mean_acceptance_upper = [ -40.4 -41.9 -45.1 -52.5 -44.7 -24.0 -16.8  -9.5  -6.5  -4.3  -3.1  -2.0 ];
log_relative_efficiency = [  -0.7  -2.7  -6.3 -14.3  -7.2  12.9  19.4  26.1  28.9  30.7  31.4  31.7 ];
log_relative_efficiency_lower = [  -0.7  -2.7  -6.3 -14.3 -23.1  -0.5  16.6  23.3  27.4  29.5  30.2  30.6 ];
log_relative_efficiency_upper = [  -0.5  -0.1   0.5   0.2  24.4  44.9  51.4  58.0  60.7  62.7  63.4  63.8 ];

% 6956 iterations
log_instantaneous_acceptance = -41.6;
nsteps_list = [      1      2      4      8     16     32     64    128    256    512   1024   2048 ];
log_mean_acceptance = [ -41.6 -43.1 -46.3 -53.7 -46.0 -25.2 -18.0 -10.5  -6.7  -4.5  -3.1  -2.1 ];
log_mean_acceptance_lower = [ -73.5 -74.2 -74.2 -71.2 -60.2 -38.3 -20.5 -12.7  -7.4  -4.7  -3.2  -2.1 ];
log_mean_acceptance_upper = [ -40.5 -42.0 -45.2 -52.6 -44.9 -24.1 -16.9  -9.6  -6.3  -4.3  -3.0  -2.0 ];
log_relative_efficiency = [  -0.7  -2.7  -6.3 -14.3  -7.2  12.9  19.4  26.2  29.3  30.8  31.5  31.9 ];
log_relative_efficiency_lower = [  -0.7  -2.7  -6.3 -14.3 -17.4   0.2  16.6  23.5  27.9  29.6  30.3  30.7 ];
log_relative_efficiency_upper = [  -0.5  -0.1   1.2   8.6  24.8  44.9  51.5  58.3  61.4  62.9  63.6  63.9 ];

% 8108 iterations
log_instantaneous_acceptance = -41.7;
nsteps_list = [      1      2      4      8     16     32     64    128    256    512   1024   2048 ];
log_mean_acceptance = [ -41.8 -43.3 -46.5 -53.8 -46.1 -25.4 -18.1 -10.7  -6.8  -4.5  -3.1  -2.1 ];
log_mean_acceptance_lower = [ -73.7 -74.3 -76.9 -71.4 -59.6 -38.4 -20.9 -12.8  -7.3  -4.7  -3.2  -2.1 ];
log_mean_acceptance_upper = [ -40.7 -42.2 -45.4 -52.7 -45.0 -24.3 -17.1  -9.7  -6.4  -4.4  -3.0  -2.0 ];
log_relative_efficiency = [  -0.7  -2.7  -6.3 -14.3  -7.2  12.9  19.4  26.2  29.4  30.9  31.7  32.0 ];
log_relative_efficiency_lower = [  -0.7  -2.7  -6.3 -14.3 -20.5   2.0  16.7  23.7  28.1  29.8  30.5  30.9 ];
log_relative_efficiency_upper = [  -0.5  -0.1   0.5   7.8  24.3  44.5  51.3  58.1  61.5  63.0  63.7  64.1 ];

% 10000 iterations
log_instantaneous_acceptance = -41.93;
nsteps_list = [      1      2      4      8     16     32     64    128    256    512   1024   2048 ];
log_mean_acceptance = [ -41.99 -43.50 -46.66 -54.06 -46.32 -25.56 -18.32 -10.76  -6.73  -4.55  -3.13  -2.07 ];
log_mean_acceptance_lower = [ -73.90 -74.50 -74.66 -71.59 -55.44 -38.46 -21.24 -12.41  -7.24  -4.71  -3.21  -2.11 ];
log_mean_acceptance_upper = [ -40.89 -42.40 -45.56 -52.96 -45.23 -24.46 -17.23  -9.91  -6.34  -4.40  -3.06  -2.03 ];
log_relative_efficiency = [  -0.75  -2.66  -6.34 -14.32  -7.22  12.87  19.44  26.31  29.65  31.14  31.87  32.23 ];
log_relative_efficiency_lower = [  -0.75  -2.66  -6.34 -14.32 -16.98   0.67  16.59  24.21  28.34  29.91  30.69  31.11 ];
log_relative_efficiency_upper = [  -0.54  -0.07   0.52   8.64  24.21  44.31  51.34  58.41  61.62  63.18  63.92  64.28 ];


%% Generate plot.

figure

scale = 1.4;

% Plot acceptance probabilities.
subplot(2,1,1);

xmin = 0.5;
xmax = nsteps_list(end) * 2;
ymin = log_instantaneous_acceptance - 1;
ymax = 0.0;

%ymin = exp(min(log_mean_acceptance) - 15);
%ymax = 1.0;

ymin = 1e-35;
ymax = 1e-0;

%plot(nsteps_list, log_mean_acceptance, 'k.');
loglog(nsteps_list, exp(log_mean_acceptance), 'k.', 'markersize', 12);

% error bars
hold on;
%loglog([nsteps_list;nsteps_list], [exp(log_mean_acceptance+2*log_mean_acceptance_error);exp(log_mean_acceptance-2*log_mean_acceptance_error)], 'k-');
loglog([nsteps_list;nsteps_list], [exp(log_mean_acceptance_upper);exp(log_mean_acceptance_lower)], 'k-');

axis([xmin xmax ymin ymax]);
hold on;
loglog([xmin xmax], exp(log_instantaneous_acceptance) * [1 1], 'r-');
xlabel('NCMC switching steps');
ylabel('acceptance probability');
set(gca, 'XTick', nsteps_list)
set(gca, 'YTick', [1e-30, 1e-20, 1e-10, 1e-0]);

% Plot efficiency boost.
subplot(2,1,2);
%plot(nsteps_list, log_relative_efficiency, 'k.');
loglog(nsteps_list, exp(log_relative_efficiency), 'k.', 'markersize', 12);
%xmin = 0;
%xmax = nsteps_list(end) * 2;
%ymin = min(log_relative_efficiency) - 1;
%ymax = max(log_relative_efficiency) + 1;

% error bars
hold on;
loglog([nsteps_list;nsteps_list], [exp(log_relative_efficiency_upper);exp(log_relative_efficiency_lower)], 'k-');

xmin = 1;
xmin = 0.5;
xmax = nsteps_list(end) * 2;
%ymin = exp(min(log_relative_efficiency) - 1);
%ymax = exp(max(log_relative_efficiency) + 10);
ymin = 1e-10;
ymax = 1e+30;

axis([xmin xmax ymin ymax]);
hold on;
loglog([xmin xmax], [1 1], 'r-');

xlabel('NCMC switching steps');
ylabel('relative efficiency gain');
set(gca, 'XTick', nsteps_list)
set(gca, 'YTick', [1e-10, 1e0, 1e+10, 1e20, 1e30]);

% Generate PDF.
filename = 'wca-dimer-efficiency.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*3.5, 'renderer', 'painters');
system(sprintf('epstopdf %s', filename));

