%% Fetch a bunch of models
mDKDP = LoadModel('Ferrell_MAPK_DKDP.txt');
mDKPP = LoadModel('Ferrell_MAPK_DKPP.txt');
mPKDP = LoadModel('Ferrell_MAPK_PKDP.txt');
mPKPP = LoadModel('Ferrell_MAPK_PKPP.txt');

m = [mDKDP; mDKPP; mPKDP; mPKPP];
nTop = numel(m);

clear mDKDP mDKPP mPKDP mPKPP

%% Simulate many models
N = 1000;

trueAll = zeros(N,1);
trueTAll = cell(N,1);
mc_pmyAll = zeros(nTop,N);
linear_pmyAll = zeros(nTop,N);
AICAll = zeros(nTop,N);
AICcAll = zeros(nTop,N);
SICAll = zeros(nTop,N);
mc_uncAll = zeros(nTop,N);
mc_weightsAll = cell(N,1);

% Warning: this loop must be parallelized on a cluster if you want this to
% finish in a reasonable amount of time
for i = 1:N
    [out{1:8}] = submitDrawMAPKTopologies(i);
    trueAll(i) = out{1};
    trueTAll{i} = out{2};
    mc_pmyAll(:,i) = out{3};
    linear_pmyAll(:,i) = out{4};
    AICAll(:,i) = out{5};
    AICcAll(:,i) = out{6};
    SICAll(:,i) = out{7};
    mc_uncAll(:,i) = out{8}.PosteriorUnc;
    mc_weightsAll{i} = out{8}.AllWeights;
end

% Note that if some iterations fail, trueAll will be zero. This is used to
% filter the results to only use the completed iterations. For various
% reasons, iterations  failed and are not included in the paper.
finished = logical(trueAll);

clear i out

%% Bootstrap estimate Monte Carlo uncertainty
mc_bootstrap_uncAll = zeros(nTop,N);
for i = 1:N
    if ~isempty(mc_weightsAll{i})
        rng(i)
        mc_bootstrap_uncAll(:,i) = bootstrap_estimate_mc_unc(mc_weightsAll{i}, [.25;.25;.25;.25], 100);
    end
end

clear i

%% Heuristic probability
lik_pmyAll = zeros(nTop,N);
aic_pmyAll = zeros(nTop,N);
sic_pmyAll = zeros(nTop,N);
for i = 1:N
    if finished(i)
        nk = vec([m.nk]);
        
        log_l = (AICAll(:,i) - 2 * nk) / -2;
        log_l = log_l - max(log_l);
        log_l = exp(log_l);
        log_l = log_l ./ sum(log_l);
        lik_pmyAll(:,i) = log_l;
        
        aic = AICAll(:,i) / -2;
        aic = aic - max(aic);
        aic = exp(aic);
        aic = aic ./ sum(aic);
        aic_pmyAll(:,i) = aic;
        
        sic = SICAll(:,i) / -2;
        sic = exp(sic);
        sic = sic ./ sum(sic);
        sic_pmyAll(:,i) = sic;
    end
end

clear i nk log_l aic sic

%% Example results
figure(1)
index = 0;
for i = [9,2,14,6]
    index = index + 1;
    subplot(2,2,index)
    bar([mc_pmyAll(:,i), linear_pmyAll(:,i), lik_pmyAll(:,i), aic_pmyAll(:,i), sic_pmyAll(:,i)])
    hold on
    errorb([1:4] - 0.305, mc_pmyAll(:,i), mc_bootstrap_uncAll(:,i), 'barwidth', .42)
    hold off
    ylim([0,1])
    xlim([0.4,4.6])
    xlabel('Topology')
    ylabel('Probability')
    if i == 9
        legend('Monte Carlo', 'Linear', 'Likelihood', 'AIC', 'SIC')
    end
end

clear i index

%% Percent correct
figure(2)

[~, mc_pmy_max] = max(mc_pmyAll,[],1);
mc_pmy_max = vec(mc_pmy_max);
mc_pmy_correct = (trueAll == mc_pmy_max);
n_mc_pmy_correct = nnz(mc_pmy_correct);

[~, linear_pmy_max] = max(linear_pmyAll,[],1);
linear_pmy_max = vec(linear_pmy_max);
linear_pmy_correct = (trueAll == linear_pmy_max);
n_linear_pmy_correct = nnz(linear_pmy_correct);

[~, lik_min] = min(AICAll,[],1);
lik_min = vec(lik_min);
lik_correct = (trueAll == lik_min);
n_lik_correct = nnz(lik_correct);

[~, AIC_min] = min(AICAll,[],1);
AIC_min = vec(AIC_min);
AIC_correct = (trueAll == AIC_min);
n_AIC_correct = nnz(AIC_correct);

[~, AICc_min] = min(AICcAll,[],1);
AICc_min = vec(AICc_min);
AICc_correct = (trueAll == AICc_min);
n_AICc_correct = nnz(AICc_correct);

[~, SIC_min] = min(SICAll,[],1);
SIC_min = vec(SIC_min);
SIC_correct = (trueAll == SIC_min);
n_SIC_correct = nnz(SIC_correct);

fraction_correct = [n_mc_pmy_correct; n_linear_pmy_correct; n_lik_correct; n_AIC_correct; n_SIC_correct;  0.25*nnz(finished)] / nnz(trueAll);

bar(fraction_correct, 'FaceColor', 'red', 'LineStyle', 'none')
hold on
errorb(fraction_correct, sqrt(nnz(finished) .* fraction_correct .* (1 - fraction_correct)) ./ nnz(finished))
hold off
ylabel('Fraction correct for most likely')
xticklabels({'MCMC', 'Linear', 'Likelihood', 'AIC', 'SIC', 'Random'})

%% JS divergence
figure(3)

mc_distances = zeros(N,1);
linear_distances = zeros(N,1);
aic_distances = zeros(N,1);
sic_distances = zeros(N,1);
lik_distances = zeros(N,1);

for i = 1:N
    if finished(i)
        mc_distances(i)     = jsdiv(mc_pmyAll(:,i), bootstrap_resample(mc_weightsAll{i}, [0.25;0.25;0.25;0.25]));
        linear_distances(i) = jsdiv(mc_pmyAll(:,i), linear_pmyAll(:,i));
        aic_distances(i)    = jsdiv(mc_pmyAll(:,i), aic_pmyAll(:,i));
        sic_distances(i)    = jsdiv(mc_pmyAll(:,i), sic_pmyAll(:,i));
        lik_distances(i)    = jsdiv(mc_pmyAll(:,i), lik_pmyAll(:,i));
    end
end

% Plot
bins = linspace(0, 1, 50);
mc_val = histc(mc_distances, bins);
linear_val = histc(linear_distances, bins);
aic_val = histc(aic_distances, bins);
sic_val = histc(sic_distances, bins);
lik_val = histc(lik_distances, bins);

mc_median = median(mc_distances(~isnan(mc_distances)));
linear_median = median(linear_distances(~isnan(linear_distances)));
lik_median = median(lik_distances(~isnan(lik_distances)));
aic_median = median(aic_distances(~isnan(aic_distances)));
sic_median = median(sic_distances(~isnan(sic_distances)));

% subplot(3,2,1)
% h = bar(bins, mc_val, 'histc');
% title('Monte Carlo')
% xlabel('JS Divergence')
% ylabel('Count')
% xlim([0 1])
% ylim([0 1000])
% set(h, 'FaceColor', 'red', 'LineStyle', 'none')
% vline(mc_median, 'k--')
% text(mc_median + 0.05, 900, sprintf('median = %0.4f', mc_median))

subplot(2,2,1)
h = bar(bins, linear_val, 'histc');
title('Linear')
xlabel('JS Divergence')
xlim([0 1])
ylim([0 1000])
set(h, 'FaceColor', 'red', 'LineStyle', 'none')
vline(linear_median, 'k--')
text(linear_median + 0.05, 900, sprintf('median = %0.3f', linear_median))

subplot(2,2,2)
h = bar(bins, lik_val, 'histc');
title('Likelihood')
xlabel('JS Divergence')
ylabel('Count')
xlim([0 1])
ylim([0 1000])
set(h, 'FaceColor', 'red', 'LineStyle', 'none')
vline(lik_median, 'k--')
text(lik_median + 0.05, 900, sprintf('median = %0.3f', lik_median))

subplot(2,2,3)
h = bar(bins, aic_val, 'histc');
title('AIC')
xlabel('JS Divergence')
xlim([0 1])
ylim([0 1000])
set(h, 'FaceColor', 'red', 'LineStyle', 'none')
vline(aic_median, 'k--')
text(aic_median + 0.05, 900, sprintf('median = %0.2f', aic_median))

subplot(2,2,4)
h = bar(bins, sic_val, 'histc');
title('SIC')
xlabel('JS Divergence')
ylabel('Count')
xlim([0 1])
ylim([0 1000])
set(h, 'FaceColor', 'red', 'LineStyle', 'none')
vline(sic_median, 'k--')
text(sic_median + 0.05, 900, sprintf('median = %0.2f', sic_median))

clear i h

%% Bias
figure(4)

mc_mean = mean(mc_pmyAll(:,finished), 2);
linear_mean = mean(linear_pmyAll(:,finished), 2);
lik_mean = mean(lik_pmyAll(:,finished), 2);
aic_mean = mean(aic_pmyAll(:,finished), 2);
sic_mean = mean(sic_pmyAll(:,finished), 2);

% Plot
subplot(3,2,1)
h = bar(mc_mean, 0.1);
title('Monte Carlo', 'Position', [2.5 0.85])
xlabel('Topology')
ylabel('Probability')
xlim([0.4 4.6])
ylim([0 1])
set(h, 'FaceColor', 'black', 'LineStyle', 'none')
hline(0.25, 'k--')

subplot(3,2,2)
h = bar(linear_mean, 0.1);
title('Linear', 'Position', [2.5 0.85])
xlabel('Topology')
ylabel('Probability')
xlim([0.4 4.6])
ylim([0 1])
set(h, 'FaceColor', 'black', 'LineStyle', 'none')
hline(0.25, 'k--')

subplot(3,2,3)
h = bar(lik_mean, 0.1);
title('Likelihood', 'Position', [2.5 0.85])
xlabel('Topology')
ylabel('Probability')
xlim([0.4 4.6])
ylim([0 1])
set(h, 'FaceColor', 'black', 'LineStyle', 'none')
hline(0.25, 'k--')

subplot(3,2,4)
h = bar(aic_mean, 0.1);
title('AIC', 'Position', [2.5 0.85])
xlabel('Topology')
ylabel('Probability')
xlim([0.4 4.6])
ylim([0 1])
set(h, 'FaceColor', 'black', 'LineStyle', 'none')
hline(0.25, 'k--')

subplot(3,2,5)
h = bar(sic_mean, 0.1);
title('SIC', 'Position', [2.5 0.85])
xlabel('Topology')
ylabel('Probability')
xlim([0.4 4.6])
ylim([0 1])
set(h, 'FaceColor', 'black', 'LineStyle', 'none')
hline(0.25, 'k--')

%% Probability comparison
figure(5)

bin = (0:0.1:1).';
%bin = logspace(-4,0,100);

subplot(3,2,1)
probcount = zeros(numel(bin),1);
truecount = zeros(numel(bin),1);
for i = 1:N
    if finished(i)
        probcount = probcount + histc(mc_pmyAll(:,i), bin);
        truecount = truecount + vec(histc(mc_pmyAll(trueAll(i),i), bin));
    end
end

% Remove junk entry from histc
probcount(end-1) = probcount(end-1) + probcount(end);
truecount(end-1) = truecount(end-1) + truecount(end);
probcount(end) = [];
truecount(end) = [];

% Plot
binprob = truecount ./ probcount;
h = bar(0.1:0.1:1, binprob);
hold on
errorb(0.1:0.1:1, binprob, sqrt(probcount .* binprob .* (1 - binprob)) ./ probcount)
hold off
set(h,'FaceColor', 'r', 'LineStyle', 'none')
title('Monte Carlo')
xlabel('Calculated probability')
ylabel('Actual probability')
xlim([0 1.1])
xtick(0.1:0.1:1)

subplot(3,2,2)
probcount = zeros(numel(bin),1);
truecount = zeros(numel(bin),1);
for i = 1:N
    if finished(i)
        probcount = probcount + histc(linear_pmyAll(:,i), bin);
        truecount = truecount + vec(histc(linear_pmyAll(trueAll(i),i), bin));
    end
end

% Remove junk entry from histc
probcount(end-1) = probcount(end-1) + probcount(end);
truecount(end-1) = truecount(end-1) + truecount(end);
probcount(end) = [];
truecount(end) = [];

% Plot
binprob = truecount ./ probcount;
h = bar(0.1:0.1:1, binprob);
hold on
errorb(0.1:0.1:1, binprob, sqrt(probcount .* binprob .* (1 - binprob)) ./ probcount)
hold off
set(h,'FaceColor', 'r', 'LineStyle', 'none')
title('Linear')
xlabel('Calculated probability')
ylabel('Actual probability')
xlim([0 1.1])
xtick(0.1:0.1:1)

subplot(3,2,3)
probcount = zeros(numel(bin),1);
truecount = zeros(numel(bin),1);
for i = 1:N
    if finished(i)
        probcount = probcount + histc(lik_pmyAll(:,i), bin);
        truecount = truecount + vec(histc(lik_pmyAll(trueAll(i),i), bin));
    end
end

% Remove junk entry from histc
probcount(end-1) = probcount(end-1) + probcount(end);
truecount(end-1) = truecount(end-1) + truecount(end);
probcount(end) = [];
truecount(end) = [];

% Plot
binprob = truecount ./ probcount;
h = bar(0.1:0.1:1, binprob);
hold on
errorb(0.1:0.1:1, binprob, sqrt(probcount .* binprob .* (1 - binprob)) ./ probcount)
hold off
set(h,'FaceColor', 'r', 'LineStyle', 'none')
title('Likelihood')
xlabel('Calculated probability')
ylabel('Actual probability')
xlim([0 1.1])
xtick(0.1:0.1:1)

subplot(3,2,4)
probcount = zeros(numel(bin),1);
truecount = zeros(numel(bin),1);
for i = 1:N
    if finished(i)
        probcount = probcount + histc(aic_pmyAll(:,i), bin);
        truecount = truecount + vec(histc(aic_pmyAll(trueAll(i),i), bin));
    end
end

% Remove junk entry from histc
probcount(end-1) = probcount(end-1) + probcount(end);
truecount(end-1) = truecount(end-1) + truecount(end);
probcount(end) = [];
truecount(end) = [];

% Plot
binprob = truecount ./ probcount;
h = bar(0.1:0.1:1, binprob);
hold on
errorb(0.1:0.1:1, binprob, sqrt(probcount .* binprob .* (1 - binprob)) ./ probcount)
hold off
set(h,'FaceColor', 'r', 'LineStyle', 'none')
title('AIC')
xlabel('Calculated probability')
ylabel('Actual probability')
xlim([0 1.1])
xtick(0.1:0.1:1)

subplot(3,2,5)
probcount = zeros(numel(bin),1);
truecount = zeros(numel(bin),1);
for i = 1:N
    if finished(i)
        probcount = probcount + histc(sic_pmyAll(:,i), bin);
        truecount = truecount + vec(histc(sic_pmyAll(trueAll(i),i), bin));
    end
end

% Remove junk entry from histc
probcount(end-1) = probcount(end-1) + probcount(end);
truecount(end-1) = truecount(end-1) + truecount(end);
probcount(end) = [];
truecount(end) = [];

% Plot
binprob = truecount ./ probcount;
h = bar(0.1:0.1:1, binprob);
hold on
errorb(0.1:0.1:1, binprob, sqrt(probcount .* binprob .* (1 - binprob)) ./ probcount)
hold off
set(h,'FaceColor', 'r', 'LineStyle', 'none')
title('SIC')
xlabel('Calculated probability')
ylabel('Actual probability')
xlim([0 1.1])
xtick(0.1:0.1:1)

clear i bin h
